AD 


Award  Number:  W81XWH-06- 1-0431 


TITLE:  Optimization  of  Tomosynthesis  Imaging  for  Improved  Mass  and 
Microcalcification  Detection  in  the  Breast 


PRINCIPAL  INVESTIGATOR:  Dan  Xia 


CONTRACTING  ORGANIZATION:  The  University  of  Chicago 

Chicago,  IL,  60637 


REPORT  DATE:  April  2009 


TYPE  OF  REPORT:  Annual  Summary 


PREPARED  FOR:  U.S.  Army  Medical  Research  and  Materiel  Command 
Fort  Detrick,  Maryland  21702-5012 


DISTRIBUTION  STATEMENT:  Approved  for  Public  Release; 

Distribution  Unlimited 


The  views,  opinions  and/or  findings  contained  in  this  report  are  those  of  the  author(s)  and 
should  not  be  construed  as  an  official  Department  of  the  Army  position,  policy  or  decision 
unless  so  designated  by  other  documentation. 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and  maintaining  the 
data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing 
this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202- 
4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently 
valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1.  REPORT  DATE 

01-04-2008 

2.  REPORT  TYPE 

Annual  Summary 

3.  DATES  COVERED 

15  Mar  2006-14  Mar  2008 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

Optimization  of  Tomosynthesis  Imaging  for  Improved  Mass  and  Microcalcification 
Detection  in  the  Breast 

5b.  GRANT  NUMBER 

W81XWH-06-1-0431 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

Dan  Xia 

5e.  TASK  NUMBER 

Email:  DANXIAtalUCHICAGO.EDU 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

The  University  of  Chicago 

Chicago,  IL,  60637 

8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

U.S.  Army  Medical  Research  and  Materiel  Command 

Fort  Detrick,  Maryland  21702-5012 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approved  for  Public  Release;  Distribution  Unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  abstract  The  goal  of  this  research  is  to  obtain  systematic  understandings  of  the  effects  of  various  physical  factors  that  are 
important  in  breast  tomosynthesis  imaging  and  to  develop  techniques  for  effectively  dealing  with  their  effects  and  for  reducing 
radiation  dose.  During  the  second  year  of  the  project  we  have  achieved  fruitful  results  based  upon  the  progress  made  in  our 
first  year  of  the  project.  Specifically,  we  have  further  investigated  the  performance  of  the  total-variation  (TV)  based  algorithm 
under  different  data  conditions  and  different  constraint  parameters.  Furthermore,  we  have  also  proposed  and  investigated  a 
new  tomosynthesis  imaging  method  with  non-planar  trajectories  for  yielding  more  data  information  with  the  same  amount  of 
imaging  dose.  We  have  also  simulated  the  scatter  in  tomosynthesis  breast  imaging  by  convolving  the  ideal  projection  data  with 
angular  dependent  scatter  kernel. 

15.  SUBJECT  TERMS 

tomosynthesis,  iterative  algorithms, 

convergence,  scanning  configuration,  physical  factors 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT 

18.  NUMBER 

OF  PAGES 

19a.  NAME  OF  RESPONSIBLE  PERSON 

USAMRMC 

a.  REPORT 

u 

b.  ABSTRACT 

u 

c.  THIS  PAGE 

u 

uu 

19 

19b.  TELEPHONE  NUMBER  (include  area 
code) 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39.18 


Table  of  Contents 


Page 

Introduction .  5 

Body .  6 

Key  Research  Accomplishments .  15 

Reportable  Outcomes .  16 

Conclusion .  17 

References .  18 


Appendices 


20 


INTRODUCTION 


Breast  tomosynthesis  is  a  tomographic  imaging  technique,  and  it  has  the  potential  advantage 
to  overcome  a  major  limitation  of  conventional  mammography  through  recovering,  to  a  large  de¬ 
gree,  the  loss  of  3D  information  about  the  breast  in  conventional,  2D  mammography  [1],  In  the 
last  several  years,  there  has  been  renewed  interest  in  developing  breast  tomosynthesis  for  detec¬ 
tion  of  breast  cancer  [2,3,4],  Although  considerable  progress  has  been  made,  improvements  to 
several  areas  of  breast  tomosynthesis  technology  are  still  needed  before  it  becomes  suitable  for 
routine  clinical  use.  In  essence,  breast  tomosynthesis  can  be  considered  as  a  dedicated  com¬ 
puted  tomography  with  limited  view  for  breast  imaging,  and  it  thus  requires  the  development  of 
special  reconstruction  algorithms  for  recovering  3D  breast  images  from  tomosynthesis  data.  In 
addition,  various  physical  factors  in  breast  tomosynthesis  can  affect  the  resulting  image  quality, 
and  the  issue  of  patient  radiation  dose  in  breast  tomosynthesis  is  of  a  concern.  The  overall  ob¬ 
jective  of  this  project  is  to  investigate  and  develop  reconstruction  algorithms  for  obtaining  breast 
images  of  practical  use,  to  investigate  and  evaluate  systematically  the  effects  of  various  physical 
factors  on  image  quality  in  breast  tomosynthesis,  and  to  use  and  evaluate  (empirical)  techniques 
for  effectively  compensating  for  the  effects  on  breast  tomosynthesis  images  and  for  possibly  re¬ 
ducing  imaging  radiation  dose  in  breast  tomosynthesis.  It  is  fully  expected  that  the  research  will 
contribute  to  the  effort  in  the  field  to  develop  and  improve  breast  tomosynthesis  for  its  clinical  use. 
This  report  summarizes  the  progress  of  this  Predoctoral  Traineeship  Award  project  made  by  the 
recipient  during  the  past  three  years. 
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BODY 


1  Training  Accomplishments 


At  the  time  of  this  report,  the  recipient  of  the  Predoctoral  Trainee- 
ship  Award  has  taken  22  out  of  22  required  courses  towards  his 
Ph.D.  degree  in  medical  physics.  The  courses  include  physics 
of  medical  imaging,  physics  of  radiation  therapy,  mathematics  for 
medical  physicists,  image  processing,  statistics,  anatomy  of  the 
body,  radiation  biology  and  teaching  assistant  training. 

2  Research  Accomplishments 

In  tomosynthesis,  data  are  acquired  only  at  a  small  number 
of  projection  views  over  a  limited  angular  range.  Therefore, 
tomosynthesis  data  are  highly  sparse  as  compared  to  data 
acquired  in  conventional  computed  tomography  (CT).  Conse¬ 
quently,  existing  analytical  algorithms  for  accurate  reconstruc¬ 
tion  of  CT  images  are  generally  not  suitable  for  yielding  useful 
images  from  tomosynthesis  data.  For  example,  the  effect  of  fil¬ 
tering  may  not  completely  be  canceled  out  due  to  the  limited 
number  of  views  and  limited  angular  range,  resulting  in  prominent  artifacts  in  reconstructed  im¬ 
ages.  In  contrast,  iterative  algorithms  can  generally  yield  images  with  less  artifacts  than  can  the 
analytic  algorithms.  In  this  project,  we  have  investigated  and  developed  a  total  variation  (TV) 
based  iterative  algorithm  for  reconstructing  accurate  images  from  incomplete  projection  data. 

The  TV  algorithm  that  is  considered  in  the  project  seeks  to  find 
the  solution  for  the  constrained  optimization  problem  [5]: 

/*  =  argmin||/j|Ty,  (1) 

subject  to  two  constraints, 

\Mf-g\<e  and  f>  0, 

where  vectors  f*  and  g  indicate  the  reconstructed  image  and  mea¬ 
sured  tomosynthesis  data,  and  matrix  M  is  the  linear  operator  rep¬ 
resenting  the  cone-beam  forward  projection  in  breast  tomosynthe¬ 
sis.  The  inequality  used  in  the  first  constraint  accounts  for  data 
inconsistency,  such  as  noise,  continuous-to-discrete  inconsistency, 
etc.  The  parameter  e  can  be  selected  for  controlling  the  impact  level 
of  potential  data  inconsistency  on  the  image  reconstruction. 


2.1  Investigation  of  reconstruction  algorithms  for  ideal 
tomosynthesis  data 

In  an  attempt  to  evaluate  the  performance  upper-bound  of  the  TV  algorithm,  we  have  first  inves¬ 
tigated  the  performance  of  the  TV  algorithm  when  applied  to  ideal  tomosynthesis  data,  which  are 
generated  from  discrete  images.  We  generated  projection  data  with  20  views  over  an  angular 
range  27r.  In  this  simulation  study,  one  has  a  discrete  matrix  system  (i.e.,  the  matrix  M  in  Eq.  (1)), 
which  allows  possible  recovery  of  the  underlying  discrete  images.  Therefore,  the  results  obtained 


Figure  2:  The  phantom  im¬ 
age  (a)  and  the  images  re¬ 
constructed  with  different  con¬ 
straint  parameter  e=1 3. 1  (b), 
e=1 1 .0  (c),  and  e=8.96  (d). 
The  display  gray  scale  is  [0.8, 
1.2], 


Figure  1:  The  Shepp-Logan 
phantom  image  (a)  and  the  im¬ 
ages  reconstructed  by  use  of  (b) 
POCS,  (c)  EM,  and  (d)  TV  algo¬ 
rithms.  The  display  gray  scale  is 
[0.9,  1.1]. 
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from  ideal  data  provide  us  the  upper  bounds  on  the  performance  of  the  TV  algorithm.  For  the  com¬ 
parison,  we  have  also  implemented  two  widely  used  existing  iterative  algorithms,  projection  onto 
convex  set  (POCS)  algorithm  [6]  and  expectation-maximization  (EM)  algorithm  [7],  The  results 
shown  in  Fig.  1  demonstrate  that,  under  ideal  situation,  the  TV  algorithm  can  yield  images  with 
superior  quality  than  POCS  and  EM  algorithms.  The  results  obtained  are  consistent  with  those 
obtained  previously  by  our  group  [5]. 

2.2  Investigation  of  image  reconstruction  inconsistent  tomosynthesis  data 

In  the  study  described  above,  we  have  investigated  the  performance  of  the  TV  algorithm  when 
applied  to  ideal  tomosynthesis  data,  which  are  generated  from  an  image  on  a  discrete  array.  How¬ 
ever,  in  realistic  tomosynthesis  imaging,  the  measured  data  are  generally  better  modeled  by  use 
of  continuous  object  functions.  Because  images  can  be  reconstructed  only  for  a  discrete  image 
from  discrete  data,  discrete  data  generated  from  a  continuous  object  function  will  contain  some 
inconsistency.  In  this  part  of  the  study,  we  carried  out  studies  on  evaluation  of  the  performance 
of  the  TV  algorithm  from  data  containing  the  continuous-to-discrete  inconsistency.  As  mentioned 
previously,  the  parameter  e  in  TV  algorithm  is  used  for  controlling  the  data  inconsistency  on  the 
image  reconstruction  and  different  constraint  parameter  e  yields  the  image  with  different  quality. 
Therefore,  selection  of  a  proper  constraint  parameter  e  is  important  to  achieve  a  good  reconstruc¬ 
tion  from  inconsistent  tomosynthesis  data. 

We  have  investigated  systematically  the  performance 
of  the  TV  algorithm  with  different  e’s.  In  one  of  the  stud¬ 
ies,  projection  data  on  20  views  over  an  angular  range 
2ir  were  generated  analytically  from  a  continuous  phan¬ 
tom  ,  as  shown  in  Fig.  2a.  Thus  the  generated  pro¬ 
jection  data  contain  the  continuous-to-discrete  inconsis¬ 
tency.  From  this  inconsistent  data  set,  images  were  then 
reconstructed  by  selecting  different  e’s.  In  Figs.  2b,  c, 
and  d,  we  show  three  images  reconstructed  by  use  of 
the  TV  algorithm  with  three  different  constraint  parame¬ 
ters  e=8.96,  11.0,  and  13.1.  The  root-mean-square  er¬ 
rors  (RMSEs)  between  the  true  image  and  three  recon¬ 
structed  images,  as  an  index  of  the  image  quality,  have 
been  calculated,  which  are  28.9,  31 .3,  and  33.2,  respec¬ 
tively.  It  can  be  observed  that,  the  RMSE  increases  as 
e  increases.  In  general,  a  large  constraint  parameter  e 
results  in  the  reconstructed  image  with  a  small  total  vari¬ 
ation  and  the  large  difference  between  the  reconstructed 
image  and  true  image. 


2.3  Improvement  of  TV  algorithm  efficiency 

The  convergence  efficiency  is  important  for  iterative  algorithms  because,  for  most  applications, 
it  is  desirable  to  obtain  sufficiently  accurate  image  in  a  given  time.  Therefore,  a  reconstruction 
algorithm  which  yields  accurate  image  within  a  10  30  iterations  is  preferable.  The  TV  algo¬ 
rithm  solves  the  constrained  minimization  problem  by  employing  adaptive  steepest  descent  POCS 
(ASD-POCS)  to  enforce  the  some  convex  constraints  on  the  image,  such  as  data  fidelity  and 
image-value  positivity,  combined  with  steepest  descent  to  reduce  TV  of  the  image.  However,  a 
large  number  of  iterations  is  required  to  reach  the  satisfactory  solution,  especially  when  the  task  is 
challenging  (detecting  a  low  contrast  object)  or  when  the  data  is  inconsistent  (noisy  or  high  scat¬ 
ter  contamination).  By  modifying  the  original  TV  algorithm,  we  have  developed  a  new  approach 


|  of  interations 

Figure  3:  The  data  distance  obtained  by 
use  of  TV  algorithm  with  different  relax¬ 
ation  parameters:  0.90  (solid  line),  0.95 
(cross),  and  0.99  (triangle).  The  image 
obtained  with  a  small  relaxation  param¬ 
eter  has  a  small  data  distance  for  each 
iteration. 
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to  improve  the  efficiency  of  TV  algorithm.  The  modification  of  the  TV  algorithm  includes  two  as¬ 
pects:  (1 )  using  ASD-POCS  with  adjusted  relaxation  parameters  to  let  data  fidelity  satisfied  within 
a  few  iterations,  and  (2)updating  the  image  along  a  new  direction,  which  is  the  combination  of  two 
gradient  directions:  data  distance  and  TV. 

We  have  investigated  systematically  the  per¬ 
formance  of  the  ASD-POCS  with  different  relax¬ 
ation  parameters.  From  the  same  projection  data 
set,  several  reconstructions  were  performed  by 
choosing  different  relaxation  parameters.  In  one 
of  these  studies,  we  generated  projection  data 
with  60  views  over  angular  range  ir.  From  this 
data  set,  we  have  reconstructed  images  from  this 
data  set  by  use  of  the  ASD-POCS  with  three 
different  relaxation  parameters,  0.90,  0.95,  and 
0.99.  The  distance  between  the  data  estimated 
from  the  reconstructed  image  and  the  generated 
data  is  used  to  indicate  the  degree  of  data  fidelity. 

Our  results  demonstrate  that  this  distance  decreases  fast  when  the  small  relaxation  parameter  is 
used,  as  shown  in  Fig.  3.  So  it  provides  a  guidance  to  adjust  the  relaxation  parameters  in  ASD- 
POCS  to  reach  the  desired  data  fidelity  constraint  within  a  small  number  of  iterations. 

The  change  of  the  image  at  each  iteration  is  designed 
to  along  a  new  direction,  which  is  formed  by  combining 
local  data  distance  gradient  direction  and  local  image  TV 
gradient  direction.  This  modified  approach  tightly  con¬ 
fines  the  updated  image  within  the  vicinity  of  the  target 
data  level  set,  and  thus  has  higher  efficiency  in  the  TV 
minimization.  We  have  modified  the  TV  algorithm  to  in¬ 
corporate  this  new  direction  for  minimizing  the  TV.  In  ex¬ 
ample  study,  we  have  generated  the  projection  data  from 
the  image  shown  in  Fig.  4a.  From  the  projection  data 
set,  the  image  was  reconstructed  by  use  of  the  proposed 
TV  algorithm,  which  is  shown  in  Fig.  4b.  For  the  com¬ 
parison,  the  image  reconstructed  by  use  of  the  original 
TV  algorithm  from  the  same  data  is  also  displayed  in  Fig. 

4c.  The  corresponding  image  TVs  are  also  displayed  in 
Fig.  5,  which  clearly  demonstrates  that  the  image  TV  ob¬ 
tained  with  the  proposed  algorithm  decreases  faster  than 
that  obtained  with  the  original  TV  algorithm.  Moreover,  the  RMSEs  of  these  two  reconstructed 
images  are  calculated,  which  are  7.25  for  new  TV  reconstruction,  8.75  for  the  original  TV  recon¬ 
struction,  respectively.  The  smaller  RMSE  obtained  with  the  new  TV  algorithm  suggests  that  at 
the  same  iteration,  the  TV  algorithm  with  the  new  designed  TV-minimization  direction  can  yield  the 
image  with  better  quality,  and,  as  much,  the  efficiency  of  the  new  TV  algorithm  has  been  improved. 

2.4  Implementation  of  TV  algorithm  to  incorporate  information  of  object  support 

The  TV  algorithm  is  a  constrained  optimization  algorithm,  which  seeks  to  find  the  image  with  mini¬ 
mum  TV  under  the  data  constraints  and  image  positivity.  Therefore,  this  algorithm  can  be  extended 
to  include  additional  constraints  such  as  object  support,  finite  lower-bound  and  upper-bound  on 
image  values.  The  incorporation  of  these  physical  constraints  to  the  reconstruction  algorithms  is 
expected  to  improve  the  image  quality.  Compared  to  structure  information  within  the  object,  infor¬ 
mation  of  the  object  support  can  be  obtained  relatively  easily  before  the  reconstruction,  therefore 


Figure  4:  The  true  image  (a)  and  images  recon¬ 
structed  by  use  of  (b)  proposed  TV  algorithm  and 
(c)  original  TV  algorithm  after  200  iterations.  No¬ 
tice  that  these  two  images  represents  the  interme¬ 
diate  images  after  certain  iterations.  The  display 
gray  scale  is  [0.9,  1.1]. 
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it  can  be  regarded  as  a  prior  information  in  our  reconstruction.  In  this  work,  we  have  developed 
an  approach  by  incorporating  the  prior  information  into  TV  reconstruction.  In  the  modified  TV  al¬ 
gorithm,  only  the  image  values  within  the  support  are  updated  while  the  image  values  outside  the 
support  is  forced  to  be  zero  during  each  iteration. 

We  have  conducted  a  study  to  investigate  the 
performance  of  the  TV  algorithm  with  the  object 
support.  We  have  generated  the  projection  data 
with  15  views  over  angular  range  27r.  The  image 
of  object  support  was  considered  as  a  prior  infor¬ 
mation,  which  is  already  known.  Two  images  are 
reconstructed  by  use  of  the  TV  with  and  without 
object  support  information,  as  shown  in  Fig.  6. 

It  can  be  observed  that  reconstruction  artifacts  in 
the  image  obtained  by  use  of  TV  algorithm  with 
object  support  information  has  been  removed. 

2.5  Investigation  of  the  effects  of  scan¬ 
ning  parameters  in  breast  tomosynthesis 

Currently,  breast  tomosynthesis  acquires  data  at 
about  12  to  20  projection  views  over  a  limited  an¬ 
gular  range  around  20  to  50  degrees.  Typically,  a 
circular  source  trajectory  is  adopted  for  collection 
of  cone-beam  projections.  It  remains,  however, 
largely  unexplored  as  to  what  the  optimal  scan¬ 
ning  configurations  are  in  tomosynthesis  imag¬ 
ing.  In  this  work,  we  have  conducted  investiga¬ 
tion  on  image  reconstruction  from  data  acquired 
at  a  small  number  of  projection  views  over  a  lim¬ 
ited  angular  range. 

In  the  study,  the  source  was  moved  along  the 
circular  trajectory.  The  projection  data  were  gen¬ 
erated  from  a  numerical  3D  breast 
phantom  in  Fig.  7a.  The  overall  shape  of  the 
breast  phantom  was  a  truncated  ellipsoid.  The 
pectoralis  muscle  was  represented  by  a  rectan¬ 
gular  slab.  The  ensemble  of  ductal  structures 
was  represented  by  a  crescent  shaped  object. 

The  mass  lesions  were  in  the  fatty  tissue,  dense 
fibroglandular  tissue,  and  the  fatty  tissue  but  with 
overlaying  dense  tissue,  respectively.  Attenu¬ 
ation  coefficients  for  breast  tissues  were  taken 
from  Ref.  [8]  and  from  the  online  table  of  NIST 
x-ray  data  for  30  keV  photons  [9]. 

We  first  generated  cone-beam  data  at  1 5,  20, 

40,  and  60  projection  views  over  27r  .  The  rea¬ 
son  for  the  angular  range  is  27r  is  that  the  study 
results  would  not  be  affected  by  the  issue  of  the 
limited  angular  range.  From  these  data  sets,  we  have  used  the  TV  algorithm  to  reconstruct  im¬ 
ages.  As  an  example,  the  images  within  2D  slices  x  —  0  and  y  =  0  reconstructed  from  15  view 
data  are  displayed  in  Fig.  7.  For  comparison,  we  have  conducted  reconstructions  by  using  POCS 


AAA 
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Figure  6:  The  true  image  (left)  and  images  recon¬ 
structed  by  use  of  TV  algorithm  with  the  objection 
support  (middle)  and  without  the  object  support 
(right).  The  display  gray  scale  is  [0.5,  1.5]. 
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Figure  7:  The  true  image  (a)  and  the  images  re¬ 
constructed  by  use  of  (b)  POCS,  (c)  EM,  and  (d) 
TV  algorithms  from  data  within  angular  ranges  2n 
for  2D  slice  x  =  0  mm  (top  row)  and  y  =  0  mm 
(bottom  row).  The  display  gray  scale  is  [0.0,  0.5]. 
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Figure  8:  The  true  image  (a)  and  images  recon¬ 
structed  from  data  within  angular  ranges  (b)  7r/4, 
(c)  7r/2,  and  (d)  n  by  use  of  TV  algorithm  for  2D 
slice  x  =  0  mm  (top  row)  and  y  =  0  mm  (bottom 
row).  The  display  gray  scale  is  [0.0,  0.5]. 
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and  EM  algorithms  as  well.  The  evaluation  indicates  that  the  TV  algorithm  is  more  accurate  than 
the  EM  and  ART  algorithms  in  reconstructing  images  from  few-view  data  and  that  images  obtained 
with  TV  algorithm  contains  less  artifacts. 

We  also  studied  image  reconstruction  from  data  acquired  over  a  limited  angular  range.  Specif¬ 
ically,  we  have  selected  angular  ranges  of  7r/4,  7t/2,  and  n  and  collected  data  at 


15,  20,  40,  and  60  views  over  each  of  these  an¬ 
gular  ranges.  From  the  collected  data,  we  re¬ 
constructed  images  by  use  of  TV,  POCS,  and 
EM  algorithms.  The  images  reconstructed  from 
15  projection  views  within  the  angular  range  7r/4, 
7t/2,  and  7r  by  use  of  TV  algorithm  are  displayed 
in  Fig.  8.  We  can  observe  that,  from  few-few 
data  over  limited  angular  range,  TV  algorithm  can 
yield  images  with  useful  information  (the  objects 
represent  the  breast  lesions  are  visible).  More¬ 
over,  these  results  suggest  that,  when  the  angu¬ 
lar  range  increases,  the  distortion  to  the  lesions  is 
suppressed,  as  shown  in  the  bottom  row  in  Fig. 8. 


2.6  Investigation  of  the  effect  of  source  trajectories  in  tomosynthesis 

x=0  y=0 

In  conventional  tomosynthesis,  the  X-ray  source 


phantom 


generally  is  moved  along  a  curve  which  is  part 
of  a  circular  trajectory.  Because  of  the  limitation 
of  the  angular  coverage,  it  is  always  challenging 
to  obtain  the  images  of  high  quality  from  con¬ 
ventional  tomosynthesis  data.  In  this  project,  we 
have  proposed  a  new  imaging  strategy  by  using 
different  source  trajectories  for  increasing  data 
information  in  tomosynthesis,  which  may  result 
in  the  improvement  of  the  reconstructed  image 
quality  without  increasing  the  imaging  radiation 
dose  [10].  In  the  investigation  of  this  new  imag¬ 
ing  approach  to  tomosynthesis  imaging,  we  have 
considered  two  new  trajectories:  one  consists  of 
two  orthogonal  arcs  and  the  other  is  a  circular 
trajectory  above  the  scanned  object,  as  shown 
in  Fig.  9.  These  two  trajectories  are  referred  to 
as  double-arc  trajectory  and  oblique-circle  trajec¬ 
tory,  respectively.  The  source  to  rotation  center 
is  7.0  cm.  A  flat-panel  detector  is  placed  per¬ 
pendicular  to  the  line  connecting  the  source  and 
the  rotation  center.  The  source-to-detector  dis¬ 
tance  is  10.0  cm.  We  have  generated  analytically 
the  data  so  that  the  data  contain  the  continuous- 
to-discrete  inconsistency.  Images  were  recon¬ 
structed  by  use  of  TV  algorithm  from  data  ac¬ 
quired  with  these  trajectories. 

In  an  example  study  of  a  series  of  investigations,  we  have  generated  projection  data  from 
the  phantom  shown  in  the  first  row  of  Fig.  10  at  30  views  uniformly  distributed  on  the  part  of  a 
circular  trajectory  used  in  breast  tomosynthesis,  whose  angular  range  is  60°.  In  the  Cartesian 


z=0 


double 

arc 


oblique 

circle 


Figure  10:  The  true  image  (first  row)  and  the 
images  reconstructed  from  projection  data  ac¬ 
quired  with  the  one-arc  trajectory  (second  row), 
the  double-arc  trajectory  (third  row),  and  oblique 
circle  trajectory  (forth  row).  The  first  to  third 
columns  represent  the  2D  slices  at  x  =  0  cm, 
y  =  0  cm,  and  z  =  0  cm,  respectively.  The  dis¬ 
play  gray  scale  is  [0,  2], 
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coordinate  system,  the  source  trajectory  can  be  written  as  r0(A)  =  ( RsinX ,  0,  RcosX),  where  A  g 
f  ].  From  the  data,  we  subsequently  reconstructed  images  by  using  the  TV  algorithm,  and  we 
displayed  in  the  second  row  of  Fig.  10  the  reconstructed  images  within  planes  specified  by  x  =  0 
cm,  y  —  0  cm,  and  z  =  0  cm. 

We  also  considered  both  double-arc  trajectory  and  oblique-circle  trajectory.  In  these  two 
scanning  configurations,  the  geometric  parameters  such  as  the  source-to-rotation-center  distance 
and  detector-to-source  distance  were  chosen  to  be  identical  to  those  in  the  above  study  for 
the  planar  circular  trajectory.  The  two  arcs  in  the  double-arc  source  trajectory  are  written  as 
ro(A)  =  (TfeinA,  0,  RcosX )  for  A  G  |]  and  ro(A)  =  (0,  i?sin(A  —  |),  Rcos(X  —  |))  for  A  G  [f ,  f], 
respectively.  The  oblique  trajectory  is  written  as  f0(A)  =  ( ifeinasinA ,  .RsinacosA,  Rcosa),  and  a  is 
a  parameter  to  describe  the  angle  between  z- axis  and  the  ray  connecting  source  and  center  of  de¬ 
tector,  which  is  30°  in  our  simulation  study.  Using  these  new  trajectories,  we  generate  analytically 
projection  data  from  the  same  phantom  over  total  30  views  of  which  are  uniformly  distributed  on 
the  trajectories.  Therefore,  the  total  imaging  radiation  doses  in  the  studies  involving  the  different 
trajectories  are  the  same.  From  the  generated  data,  we  reconstructed  images  by  use  of  the  TV 
algorithm.  In  the  third  and  forth  rows  of  Fig.  10,  we  show  the  reconstructed  images  for  these  two 
trajectories  within  planes  specified  by  x  =  0  cm,  y  =  0  cm,  and  z  =  0.5  cm,  respectively.  Com¬ 
parison  of  the  results  in  Fig.  10  suggests  that  with  the  same  number  of  views  (or,  equivalently, 
the  same  amount  of  imaging  dose),  data  acquired  with  our  proposed  trajectories  tomosynthesis 
may  contain  more  information  than  that  acquired  with  the  conventional  tomosynthesis,  leading  to 
images  with  improved  quality. 

2.7  Investigation  of  the  physical  factors  in  tomosynthesis  imaging 

In  realistic  tomosynthesis  imaging,  a  number  of  physical  factors  can  significantly  affect  image 
quality,  thereby  lowering  the  detection  rate  of  mass  and  microcalcification  in  breast.  In  this  project, 
we  have  conducted  investigation  of  the  impact  of  the  major  physical  factors,  including  data  noise, 
scatter,  non-linear  partial  volume,  beam-hardening,  and  non-isotropic  image  spatial  resolution. 
We  summarize  briefly  our  investigations  below. 

2.7.1  Noise  properties 

For  tomosynthesis  breast  imaging  that  aims  to  have  comparable 
radiation  dosage  with  the  standard  mammography,  the  resulting 
data  will  contain  substantial  noise  because  the  radiation  expo¬ 
sure  in  one,  or  two,  view  in  the  standard  mammography  will  be 
divided  among  many  more  projection  views  in  tomosynthesis.  In 
general,  images  generated  from  these  tomosynthesis  data  will 
be  noisier,  leading  to  decreased  detectability  of  microcalcifica¬ 
tion  and  low-contrast  structures.  The  noise  characteristics  of  the 
images,  however,  depend  significantly  on  the  reconstruction.  A 
systematic  noise  analysis  in  the  context  of  the  tomosynthesis 
imaging  is  required. 

In  an  attempt  to  separate  the  issues  of  few-view  and  limited 
angular  range  from  data  noise,  we  have  first  performed  a  thorough  noise  study  for  image  recon¬ 
struction  from  projection  data  collected  at  a  large  number  of  views  over  2ir  [11],  The  result  of  this 
study  provides  a  theoretical  guidance  to  the  investigation  of  the  noise  properties  in  tomosynthe¬ 
sis  in  which  data  are  acquired  at  a  small  number  of  views  over  a  limited  angular  range.  In  this 
study,  we  specifically  focus  on  investigating  how  data  noise  are  propagated  into  the  reconstructed 
images.  Both  theoretical  and  numerical  analyses  were  carried  out,  and  the  results  indicate  that 
variances  of  reconstructed  images  are  spatially  varying  and  that  the  levels  of  variances  in  different 


Figure  1 1 :  Empirical  variance  im¬ 
ages  obtained  by  use  of  ana¬ 
lytical  algorithm  from  fan-beam 
data  containing  Gaussian  (left) 
and  Poisson  (right)  noise. 
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regions  are  not  affecting  significantly  each  other,  as  shown  in  Fig.  1 1 .  Because  the  study  was 
based  upon  a  theoretical  result,  it  was  computationally  possible  to  accomplish  this  initial  study. 
On  the  other  hand,  noise  studies  involve  the  iterative  algorithms  such  as  the  POCS,  EM,  and  TV 
algorithms  are  much  more  demanding  computationally. 

2.7.2  X-ray  scatter  effect 

For  most  existing  reconstruction  algorithms,  the  projection  data 
are  assumed  to  include  only  the  x-ray  transmission  on  a  straight 
line  connecting  the  source  and  detector  element.  However,  a 
fraction  of  scattered  x-ray  intensity  is  also  involved  in  realistic 
measurements.  The  inconsistencies  introduced  by  the  presence 
of  scattered  radiations  in  the  projection  data  can  lead  to  a  de¬ 
crease  of  low-contrast  detectability,  cupping  artifacts  and  streak 
artifacts  between  dense  objects.  The  effect  of  X-ray  scatter 
poses  a  challenging  problem  in  image  reconstruction  because 
the  precise  functional  form  of  the  distribution  of  the  scattered 
radiations  depends  on  the  subject  being  scanned  [12]. 

Taking  advantage  of  the  fact  that  the  scattered  X-ray  inten¬ 
sity  does  not  display  significant  high-resolution  structures,  we 
have  proposed  a  convolution  method  in  our  study  of  scatter  ef¬ 
fect.  Specifically,  we  have  employed  an  existing  scatter  kernel  to  generate  scatter  components  in 
tomosynthesis  projection  data  [13].  In  this  model,  critical  physical  parameters  such  as  chest  wall 
to  nipple  distance,  breast  thickness,  and  breast  glandularity  have  been  taken  into  account.  We 
have  generated  and  evaluated  scatter  components  for  projection  angles:  0°,  6°,  12°,  18°,  24°, 
and  30°,  which  are  typical  angles  in  breast  tomosynthesis.  Based  upon  the  convolution  kernel  and 
scatter  components,  we  determined  the  scatter  kernel  for  any  projection  angle  by  use  of  interpo¬ 
lation/extrapolation  methods.  Once  we  obtain  these  kernels  for  any  angles,  we  then  estimated  the 
the  scatter-less  projection  data  by  deconvolving  with  the  scatter  kernel. 

In  an  attempt  to  test  the  algorithm,  we  have  applied  this  method  to  image  reconstruction  from 
data  acquired  with  a  real  breast  CT  scanner,  in  which  the  source  and  the  flat  panel  can  be  rotated 
around  the  patient’s  breast.  For  each  scan,  projection  data  were  collected  at  typical  500  views 
uniformly  distributed  over  360°.  From  the  full  data  set,  we  selected  a  subset  of  data  at  20  views 
uniformly  distributed  over  72°  for  simulating  tomosynthesis  data.  For  each  projection,  the  scatter¬ 
less  projection  were  estimated  from  the  original  data  by  using  the  deconvolution  method.  From 
the  corrected  projection  data  we  reconstructed  the  breast  image  by  use  of  TV  algorithm,  the  result 
is  shown  in  Fig.  12a.  For  comparison,  the  image  reconstruction  from  projection  data  with  scatter 
was  also  performed,  which  is  shown  in  Fig.  12b.  It  can  be  observed  that  the  contrast  in  Fig.  12a 
is  higher  than  that  in  12b,  which  suggests  that  the  compensation  of  scatter  effect  can  improve  the 
image  quality  in  tomosynthesis  imaging. 

2.7.3  Non-linear  partial  volume  effect 

The  ideal  imaging  model  in  tomosynthesis  imaging  assumes  a  point-like  x-ray  source  and  infinitely 
small  detector  elements.  However,  in  realistic  situation  the  source  and  detector  elements  both 
have  finite  sizes.  As  a  result,  the  measured  data  reflect  the  averaged  projections,  which  are 
equal  to  the  negative  exponentials  of  the  path  integrals,  over  the  areas  of  the  source  spot  and/or 
the  detector  bin.  Obviously,  this  process  for  generating  the  averaged  line  integrals  violates  the 
ideal  imaging  model.  In  CT  images,  this  non-linear  averaging  can  result  in  conspicuous  artifacts, 
such  as  streak  artifacts  that  occur  tangential  to  the  structures  with  sharp  attenuation  transitions 


Figure  12:  Images  reconstructed 
from  data  with  scatter  correction 
(left)  and  without  scatter  correc¬ 
tion  (right).  The  display  gray 
scale  is  [0.18,  0.29]. 
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[14,  15].  Therefore  in  this  study,  we  have  investigated  the  non-linear  partial  volume  effect  in  the 
tomosynthesis  imaging,  especially  for  the  detection  of  breast  microcalcifications. 

In  our  simulation,  the  projection  data  were  generated 
from  a  numerical  3D  breast  phantom  with  some  microcalci¬ 
fications,  as  shown  in  Fig.  13a.  The  smallest  white  points 
in  the  image  represent  the  microcalcifications.  We  have 
generated  tomosynthesis  data  at  30  projection  views  uni¬ 
formly  distributed  over  the  angular  range  7r/4.  For  each 
project  view,  we  have  sub-divided  the  x-ray  source  with  finite 
size  0.75  x  0.75  mm2  and  detector  elements  with  finite  size 
0.75  x  0.75  mm2  into  5x5  sub-elements  and  computed  the 
integrals  over  the  625  lines  joining  the  sub-elements  of  the 
source  and  the  sub-elements  of  the  detector  element.  The 
projection  data  for  each  element  were  then  obtained  by  tak¬ 
ing  the  logarithm  of  the  average  of  the  negative  exponentials 
of  these  line  integrals.  From  the  generated  tomosynthesis 
data,  the  image  was  reconstructed  by  use  of  TV  algorithm, 
which  is  displayed  in  Fig.  13b.  Compared  to  the  image  re¬ 
constructed  from  the  data  without  modeling  non-linear  par¬ 
tial  volume  effect,  as  shown  in  Fig.  13c,  there  is  no  obvious  additional  artifacts  in  the  image 
obtained  from  data  with  non-linear  partial  volume  effect.  The  result  suggests  that  in  tomosynthe¬ 
sis,  the  artifacts  caused  by  non-linear  partial  volume  effect  is  overwhelmed  by  the  artifacts  caused 
by  other  effects,  such  as  limited  scanning  angular  range. 

2.7.4  Beam-hardening  effect 

Beam  hardening  is  an  unavoidable  physical  factor  that  needs 
to  be  accounted  in  x-ray  systems.  In  most  existing  work  on 
algorithm  development  for  image  reconstruction,  the  projec¬ 
tion  data  are  assumed  to  be  the  transmission  of  monochro¬ 
matic  x-rays  through  the  subject.  A  realistic  x-ray  source, 
however,  has  a  broad  spectrum  and  is  not  monochromatic. 

Because  the  x-ray  attenuation  of  an  object  depends  on  the 
photon  energy,  with  lower-energy  photons  being  attenuated 
more  strongly,  the  mean  energy  of  the  x-ray  beam  will  in¬ 
crease  as  it  penetrates  through,  i.e.,  at  exit  the  beam  be¬ 
comes  harder  than  original.  Consequently,  the  monochro¬ 
matic  imaging  model  is  incorrect  and  this  effect,  without  com¬ 
pensation,  can  cause  artifacts  such  as  cupping,  streaking, 
and  an  overall  shift  of  the  CT  number  in  the  reconstructed 
images  [16, 17]. 

We  have  investigated  the  beam  hardening  effect  by  con¬ 
sidering  a  realistic  x-ray  spectrum  in  our  simulations.  In  this 
a  28  kVp  setting  with  a  molybdenum  anode.  Based  upon  the  spectrum  of  this  source,  we  obtained 
the  relative  photon  fluence  at  energy  E  =  5,  10,  15,  20  25,  30  keV,  respectively.  The  attenuation 
coefficient  within  a  numerical  breast  phantom  (Fig.  4a)  is  designed  as  a  function  of  the  x-ray  en¬ 
ergy,  which  were  taken  from  Ref.  [8].  Based  upon  the  known  relative  photon  fluence,  we  have 
accordingly  weighted  the  exponentials  of  line  integrals  for  each  individual  energy,  then  took  the 
logarithm  of  the  sum  of  the  weighted  exponentials  to  form  the  final  projection  data  for  the  poly¬ 
chromatic  source.  With  this  method,  we  generate  the  tomosynthesis  data  by  using  the  scanning 
configuration  described  in  Sec.  2.7.3.  The  image  was  reconstructed  by  use  of  the  TV  algorithm, 
which  is  shown  in  Fig.  14c.  For  the  comparison  the  image  reconstruction  for  monochromatic  x-ray 


Figure  14:  Breast  phantom  image  for 
30keV  energy  (a)  and  the  images  re¬ 
constructed  from  data  acquired  with 
30keV  monochromatic  x-ray  source 

(b)  and  realistic  x-ray  source  data 

(c) .  The  display  gray  scales  for  fig¬ 

ures  (a)  and  (b)  are  [0.1,  0.5],  while 
the  display  gray  scale  for  figure  (c)  is 
[0.2, 1.0]  because  of  the  CT  number 
shift. _ 

study,  the  x-ray  source  was  modeled 


Figure  13:  Breast  phantom  image 
with  some  microcalcifications  (a)  and 
images  reconstructed  from  data  ac¬ 
quired  with  modeling  non-linear  par¬ 
tial  volume  effect  (b)  and  data  ac¬ 
quired  without  modeling  non-linear 
partial  volume  effect  (c).  The  display 
gray  scale  is  [0.0,  0.5]. 
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was  also  performed,  the  result  is  displayed  in  Fig.  14b.  One  can  find  the  severe  artifacts,  such 
as  cupping  and  CT  number  shift,  in  the  image  reconstructed  from  data  simulated  by  using  poly¬ 
chromatic  x-ray.  Therefore,  the  correction  of  beam-hardening  in  tomosynthesis  imaging  may  be 
necessary,  which  is  now  under  investigation. 


2.7.5  Non-isotropic  resolution 

In  current  breast  tomosynthesis,  the  spatial  resolution  within  a 
transverse  plane  is  much  finer  than  that  along  the  longitudinal 
direction,  and  image  representation  with  non-isotropic  image 
size  is  used  in  iterative  algorithms  for  reducing  computational 
time.  Such  an  image  representation  can  lead  to  significant  arti¬ 
facts  in  iterative  reconstruction.  We  have  been  investigating  the 
effect  of  non-isotropic  image  representation  on  iterative  recon¬ 
struction  accuracy  of  breast  tomosynthesis  images  [18].  We 
have  reconstructed  images  by  use  of  POCS,  EM,  and  TV  al¬ 
gorithms  for  image  representations  with  different  ratios  of  the 
in-plane  and  longitudinal  pixel  size,  which  are  shown  in  Fig.  1 5. 

Our  results  demonstrate  that  non-isotropic  image  representa¬ 
tion  can  lead  to  image  artifacts  in  reconstructed  images.  The 
appearance  and  severity  of  the  artifacts  depend  not  only  upon 
the  ratio  between  the  in-plane  and  longitudinal  resolution  but 
also  upon  the  iterative  algorithms.  The  TV  algorithm  seems  to  be  less  susceptible  to  the  effect 
than  the  POCS  and  EM  algorithms. 


Phantom 


POCS 


EM 


TV 


Figure  15:  The  true  image  (a) 
and  the  images  reconstructed  with 
isotropic  pixel  size  (left)  and  non¬ 
isotropic  spatial  pixel  size  (right). 
The  display  gray  scale  is  [0.0,  0.5]. 
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KEY  RESEARCH  ACCOMPLISHMENTS 


•  We  have  modified  the  TV  algorithm  tailored  to  breast  tomosynthesis. 

•  We  have  implemented  the  existing  EM  and  POCS  algorithms  for  comparison  with  our  TV 
algorithm. 

•  We  have  determined  the  upper  bound  on  the  performance  of  TV  algorithm  for  image  recon¬ 
struction  from  discrete  tomosynthesis  data. 

•  We  have  implemented  a  projection  program  to  compute  both  discrete-to-discrete  projection 
and  continuous-to-discrete  projection  data. 

•  We  have  conducted  an  investigation  on  the  effect  of  different  constraint  parameter  e  for  the 
data  containing  continuous-to-discrete  inconsistency. 

•  We  have  modified  TV  algorithm  for  improvement  of  efficiency  by  using  different  relaxation 
parameter. 

•  We  have  proposed  an  approach  for  efficiently  minimizing  TV  in  the  reconstruction. 

•  We  have  implemented  the  TV  algorithm  to  incorporate  the  object  support  information  to 
improve  tomosynthesis  reconstruction. 

•  We  have  investigated  the  effects  of  scanning  parameters  in  breast  tomosynthesis  imaging 
by  using  different  numbers  of  views  over  different  angular  ranges. 

•  We  have  conducted  a  preliminary  investigation  of  the  different  source  trajectories  in  tomosyn¬ 
thesis  breast  imaging  for  improving  image  quality. 

•  We  have  performed  a  preliminary  investigation  of  the  data  noise. 

•  We  have  generated  the  scatter  component  in  breast  tomosynthesis  by  involving  the  scatter- 
free  projection  data  with  the  scatter  kernel. 

•  We  have  applied  a  deconvolution  method  to  real  breast  data  to  improve  the  image  recon¬ 
struction. 

•  We  have  performed  a  preliminary  investigation  on  the  effect  of  the  non-linear  partial  volume 
in  breast  tomosynthesis. 

•  We  have  conducted  a  preliminary  investigation  on  the  effect  of  beam-hardening  in  breast 
tomosynthesis. 

•  We  have  performed  a  preliminary  investigation  on  the  effect  of  the  non-isotropic  spatial  res¬ 
olution  in  breast  tomosynthesis. 
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CONCLUSIONS 


The  recipient  of  the  Predoctoral  Traineeship  Award  has  finished  the  required  courses  towards 
his  Ph.D.  degree.  These  trainings  have  proved  useful  for  the  recipient  to  achieve  the  proposed 
research  goals. 

During  the  research,  we  have  implemented  and  investigated  the  TV  algorithm  for  image  re¬ 
construction  in  breast  tomosynthesis.  We  have  conducted  numerical  studies  to  investigate  the 
image  reconstruction  from  ideal  projection  data  by  use  of  the  TV  algorithm  and  the  existing  EM 
and  TV  algorithms,  which  give  the  upper  bound  on  the  performance  of  these  reconstruction  algo¬ 
rithms.  Our  results  suggest  that,  in  general,  the  TV  algorithm  provide  tomosynthesis  images  with 
quality  higher  than  EM  and  TV  algorithms  in  terms  of  the  metrics  such  as  root-mean-square  er¬ 
ror.  We  have  performed  the  numerical  studies  to  investigate  the  image  reconstruction  for  the  data 
containing  continuous-to-discrete  inconsistency,  which  happens  in  real  scans.  In  this  study,  the 
images  were  reconstructed  by  use  of  the  TV  algorithm  with  different  constraint  parameter  e.  The 
relationship  between  the  image  quality,  in  terms  of  the  metrics  such  as  RMSE,  and  the  constraint 
parameter  e  has  been  identified.  Base  upon  the  original  TV  algorithm,  we  have  further  modified 
TV  algorithm  to  improve  the  performance  of  image  reconstruction,  including  improving  the  TV  al¬ 
gorithm  efficiency  and  incorporating  the  information  of  object  support.  The  TV  reconstruction  with 
information  of  object  support  yields  the  better  image  than  does  the  original  TV. 

In  this  project,  we  have  also  investigation  of  different  scanning  configurations  in  breast  to¬ 
mosynthesis.  The  tomosynthesis  data  were  collected  from  different  angular  ranges  at  different 
numbers  of  views.  Our  results  suggest  that  when  the  angular  range  decreases,  image  quality 
obtained  with  these  algorithms  decreases.  Moreover,  we  have  proposed  a  new  imaging  strat¬ 
egy  by  using  trajectories  for  increasing  data  information  in  tomosynthesis.  The  results  obtained 
with  the  new  scanning  trajectories  demonstrate  that  with  the  same  number  of  views  (or,  equiva¬ 
lently,  the  same  amount  of  imaging  dose),  data  acquired  with  the  proposed  non-planar  trajectory 
tomosynthesis  may  contain  more  information  than  that  acquired  with  the  conventional  tomosyn¬ 
thesis,  leading  to  images  with  improved  quality. 

Furthermore,  we  have  also  performed  the  studies  to  model  some  critical  physical  factors  in 
breast  tomosynthesis,  including  data  noise,  scatter,  non-linear  partial  volume,  beam  hardening 
and  non-isotropic  spatial  resolution  on  tomosynthesis  images.  The  effect  of  these  physical  factors 
has  been  investigated.  Our  results  demonstrate  that  noise,  scatter,  beam  hardening  can  affect 
the  image  quality  obtained  from  the  tomosynthesis  data,  which  should  be  compensated  during  the 
reconstruction  process,  while  the  artifacts  caused  by  the  non-linear  partial  volume  effect  in  breast 
tomosynthesis  is  negligible.  Moreover,  a  deconvolution  method  for  scatter  correction  has  been 
applied  for  real  breast  data,  and  has  been  demonstrated  to  be  helpful  for  improving  the  image 
quality. 

Overall,  we  have  achieved  the  goals  for  this  project  and  laid  the  foundation  for  our  continuation 
of  breast  cancer  research.  The  future  research  will  continue  on  the  evaluation  and  correction  of 
physical  factors  such  as  beam-hardening.  Moreover,  the  developed  TV  algorithm  will  be  applied 
to  real  data  of  breasts  for  additional  evaluation  of  the  scanning  configurations  and  reconstruction 
algorithms.  For  this,  real-patient  data  base  will  be  needed  for  testing  and  assessing  the  recon¬ 
struction  algorithms  and  imaging  parameters. 
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Abstract — Recently,  there  has  been  much  progress  in  algorithm 
development  for  image  reconstruction  in  cone-beam  computed  to¬ 
mography  (CT).  Current  algorithms,  including  the  chord-based  al¬ 
gorithms,  now  accept  minimal  data  sets  for  obtaining  images  on 
volume  regions-of-interest  (ROIs)  thereby  potentially  allowing  for 
reduction  of  X-ray  dose  in  diagnostic  CT.  As  these  developments 
are  relatively  new,  little  effort  has  been  directed  at  investigating 
the  response  of  the  resulting  algorithm  implementations  to  physical 
factors  such  as  data  noise.  In  this  paper,  we  perform  an  investiga¬ 
tion  on  the  noise  properties  of  ROI  images  reconstructed  by  using 
chord-based  algorithms  for  different  scanning  configurations.  We 
find  that,  for  the  cases  under  study,  the  chord-based  algorithms 
yield  images  with  comparable  quality.  Additionally,  it  is  observed 
that,  in  many  situations,  large  data  sets  contain  extraneous  data 
that  may  not  reduce  the  ROI-image  variances. 

Index  Terms — Chord,  computed  tomography  (CT),  cone-beam 
CT,  noise,  reconstruction. 


I.  Introduction 

IN  recent  years,  exact  algorithms  have  been  developed  for 
reconstructing  images  [1]  and  for  reconstructing  images  on 
“7r-lines”  [2]-[4]  from  helical  cone-beam  data.  Since  2005,  pa¬ 
pers  have  being  published  on  algorithm  development  for  re¬ 
constructing  images  on  chords  for  general  trajectories  [5]— [8] . 
Some  of  these  algorithms  can  reconstruct  images  within  3-D 
regions  of  interest  (ROIs)  from  cone-beam  data  containing  both 
longitudinal  and  transverse  truncations.  The  introduction  of  the 
M- line  concept  and  reconstruction  [5],  [9]  provides  additional 
flexibility  for  covering  volume  ROIs. 

As  these  algorithm  developments  are  relatively  recent,  little 
effort  has  been  directed  at  investigating  their  noise  properties. 
With  the  algorithm  development  for  ROI-image  reconstruction, 
it  has  been  tacitly  assumed  that  the  reduction  in  necessary  scan¬ 
ning  angle  and  in  projection  data  may  lead  to  ROI  images  from 
less  radiation  exposure.  This  conclusion  may,  however,  depend 
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on  the  noise  properties  of  reconstruction  algorithms.  If  ROI  re¬ 
construction  from  the  minimal  (or  reduced)  data  set  leads  to 
noisier  ROI  images  than  reconstruction  of  the  same  ROI  from  a 
larger  data  set,  it  may  be  necessary  to  increase  the  X-ray  source 
intensity  for  the  ROI  data  set  to  attain  the  same  image  quality 
as  those  reconstructed  from  larger  data  sets.  Such  an  increase 
can  offset  the  fact  that  reduced  or  minimum  projection  data  are 
needed  for  ROI  reconstruction. 

The  focus  of  this  paper  is  to  investigate  the  noise  proper¬ 
ties  of  image  reconstruction  from  minimal  data  set  and  large 
data  set  by  use  of  chord-based  algorithms.  We  demonstrate 
that  the  minimal  data  set  can  indeed  lead  to  actual  reduc¬ 
tion  of  radiation  exposure  for  attaining  comparable  image 
quality,  defined  in  terms  of  image  variance,  as  that  obtained 
with  a  larger  data  set.  In  Section  II,  we  briefly  summarize 
the  chord-based  reconstruction  algorithms:  backprojection  fil¬ 
tration  (BPF)  [2],  [6],  minimum  data  filtered  backprojection 
(MDFBP)  [4],  [6],  and  filtered  backprojection  (FBP)  [6],  [10] 
algorithms.  In  Section  III,  we  perform  analysis  and  empir¬ 
ical  studies  on  noise  properties  of  images  reconstructed  from 
parallel-beam,  fan-beam,  and  cone-beam  data.  Finally,  a  dis¬ 
cussion  is  given  in  Section  IV. 

II.  Chord-Based  Reconstruction  Algorithms 

We  consider  a  continuous  source  trajectory  specified  by 
r0(s)  =  (x(s),  y(s),  z(<s)),  where  x(s),  y(s ),  and  z(«s)  denote 
the  x— ,  y— ,  and  components  of  fo(s)  in  the  fixed-coordi¬ 
nate  system,  and  s  is  a  curve  parameter  indicating  the  position 
of  the  X-ray  source  on  the  trajectory.  The  projection  data  of  the 
object  function  f(f)  can  be  mathematically  expressed  as 

oo 

D  (r0(s),/3)  =  J  dtf  (ro(s)  +tjAj  (1) 

o 

where  the  unit  vector  /3  denotes  the  direction  of  a  specific  X-ray 
passing  through  the  point  r.  We  also  introduce  two  additional 
coordinate  systems  {u^v^w}  and  {uc %$Vd}  to  describe  the 
geometry  in  a  general  scan.  They  are  fixed  on  the  rotating 
source  point  and  the  cone-beam  projection  of  the  source  point, 
respectively,  which  are  referred  to  as  the  rotation-coordinate 
and  detector-coordinate  systems.  Let  eu(s),  ev(s),  and  ew(s) 
denote  the  orthogonal  unit  vectors  of  the  rotation-coordinate 
system.  The  rotation-coordinate  system  can  be  chosen  such  that 
eu(s)  and  ew(s)  are  within  the  x  —  y  plane  and  ev(«s)  is  parallel 
to  the  z  axis.  One  can  also  choose  the  “well  oriented”  coor¬ 
dinate  system  as  the  rotation-coordinate  system  [5],  in  which 
unit  vector  eu(s)  is  parallel  to  and  unit  vectors  e„(s)  and  e^s) 
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are  perpendicular  to  the  direction  of  ( dr0(s)/ds ).  We  assume 
that  a  detector  plane  is  placed  at  a  distance  S  from  the  source 
point  and  orients  along  ew(s).  The  detector-coordinate  system 
{ud,vd}  is  the  cone-beam  projection  of  the  2-D  coordinate 
system  {u,v}  onto  the  detector  plane,  and  the  ud  and  vd  axis 
are  along  eu(s)  and  6„(.sj,  respectively.  In  this  situation,  we 
also  use  P(ud,  vd,  s)  to  denote  the  cone-beam  projection,  thus 
D(r0(s),/3)  =  P(ud,vd,s),  when 


specified  by  sa  and  sb.  The  filtered  image  f(xc ,  sa,  sb )  is  given 

by 

f(Xc,  Sa,  Sb )  =  ^  J  «.  «a,  Sb)  (5) 

R 

where  the  backprojection  image  on  the  chord  can  be  written  as 


/?  = 


A{ud,vd) 
and  A{ud, vd)  =  yf^+^+S2. 


f"v/e„(-s!  +  vdev(s )  -  .‘>n„.(.s)] 


(2) 


In  a  2-D  case,  it  can  be  observed  that  vd  =  0.  For  notational 
convenience,  we  use  A{ud)  and  P{ud,  s)  to  denote  A(ud,  0)  and 
P(ud,0,s),  respectively. 

A  chord  is  a  line  segment  connecting  two  points  f0(sa )  and 
r0(sb)  on  the  trajectory.  Any  point  r  on  the  chord  can  be  ex¬ 
pressed  as 


IcSKi] 


(3) 


g(x'c,sa,sb )  =  nc  (x'c)  /  ds 


sgn(-/3  •  e^) 

I P  ~  To(s)|2 


x<(  -dr°js^  ■  (3P{ud,vd,s) 


+ 


ds 
dr0(s) 


ds  uy  '  ’  S(s)  ds 
8P(ud,vd,s) 


x  A(ud,vd) 


dud 


ds  1  ds 

9P(ud,  vd,  s) 


x  A{ud,  vd ) 


dvd 


(6) 


where  ec  =  (f0(st)-fo(sa))/|fo(sfc)-fo(sa)|  denotes  the  di¬ 
rection  of  the  chord,  and  l  =  (l/2)|f0(st)  -  fo(sa)|  is  one  half 
of  the  chord  length.  For  a  helical  trajectory,  the  curve  parameter 
s  is  linearly  related  to  the  rotation  angle  A,  and  in  the  current 
work,  we  selects  =  A.  When  sa  and  sb  are  within  one  turn,  the 
chord  becomes  the  conventional  7r-line  segment  [2],  [11],  [12], 
The  intersection  between  a  chord  and  the  object  is  referred  to  as 
a  support  segment.  Let  xc\  and  xc2  represent  the  end  points  of 
a  support  segment.  Because  the  trajectory  under  consideration 
never  intersects  the  object,  we  have  [xci,xc2]  c  [- 1,  l ].  There¬ 
fore,  one  can  use  ( xc,sa,sb )  and  fc(xc,sa,sb )  to  denotea  point 
and  the  corresponding  image  on  the  chord.  We  have  previously 
developed  three  algorithms,  which  are  referred  to  as  the  BPF 
[2],  [6],  [10],  M  DFBP  [4],  [6],  and  FBP  [6],  [10]  algorithms, 
respectively,  for  exact  image  reconstruction  on  a  chord  of  a  gen¬ 
eral  trajectory. 

A.  BPF  Algorithm 

The  BPF  algorithm  [2],  [6]  reconstructs  the  image  on  a  chord 
specified  by  sa  and  sh  as 


and  the  rect  function  Uc(x'c)  =  1  if  x'c  e  [xA,xB]  and  zero 
otherwise.  It  can  be  observed  in  (4)  that  the  chord  image  can  be 
obtained  exactly  from  knowledge  of  the  backprojection  image 
g(x'c,  sa,  sb )  for  x'c  e  [xA,  xb],  which  we  refer  to  as  the  recon¬ 
struction  segment  because  it  determines  the  actual  reconstruc¬ 
tion  interval  on  the  chord.  In  particular,  because  the  reconstruc¬ 
tion  segment  [xA,xB\  can  be  chosen  as  small  as  the  support 
segment  [: xcl,xc2 ],  the  chord  image  can  be  reconstructed  from 
knowledge  of  g(x'c,sa,sb)  only  on  the  support  segment.  This 
interesting  property  of  the  H  ilbert  transform  forms  the  basis  for 
exact  image  reconstruction  on  a  chord  from  projections  con¬ 
taining  longitudinal  or  transverse  truncations  [13], 

B.  MDFBP  Algorithm 

The  BPF  algorithm  reconstructs  the  chord  image  by  per¬ 
forming  a  1-D  filtration  [i.e.,  the  integration  over  x’c  in  (4)]  of 
the  backprojection  image  [i.e.,  the  integration  over  s  in  (6)]. 
On  the  other  hand,  the  MDFBP  algorithm  [4],  [6]  reconstructs 
the  chord  image  by  performing  a  1-D  data  filtration  (i.e.,  the 
integration  over  u'c)  prior  to  their  backprojection  (i.e.,  the 
integration  over  s)  onto  the  chord 


/cC'l'C,  ^a,  £&)  —  /('1'c,  ^a,  £&)  “P 


1  P{tldo,Vdo,  So) 

2ir  b(xc ) 


s/JT—  ®b)(^  —  Xa)  \f(l  +  £a)((  +  Fb) 


l  —  xc 


+ 


l  +  xc 


(4) 


where  xc  e  [xA,xB\,  and  parameters^  and^s  are  two  points 
on  the  chord  satisfying  [xcl,xc2]  c  [xA,xB]  c  [-/,/].  The 
function  b(xc)  is  defined  as  b(xc)  =  -J{x^  -  a;c)(®c  -  ~xa), 
and  P(udo,vdo,sa )  denotes  the  projection  along  the  chord  22 


f r ( Sa,  Sb ) 


49 b 


27r2  b(x , 

Sa  R 

1  P(UdO,VdO,Sa) 

27 r  b(xc ) 

aJ(1-  xB)(l  ~  xA)  y/(l  +  xA)(l  +  xB) 


Pu 


+ 


l  —  xc 


+ 


/  xc 


X 


(7) 
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where  the  modified  data  function  is  given  by 


Pu  =  n  c 


b(x'c)  sgn(-p-ew 

w2  (1  -  O  +  «>i<  | _  ro(s)|2 

d  *fiP{u<l,Vd,s) 

\dr0(s)  „  ud  rfr0(s) 


^(f)  .a  /'o',  ,  ^o(s) 

ds  u  S(s)  ds 


x  A(ud,Vd ) 
.  [dfoO) 


dPjUd,  V<U  g) 
dud 


"oOO  a  /  x  ,  Vd  dfojs) 
ds  v  S(s)  ds 


X  A(ud,  Vd) 


dP(ud,vd,s) 


W\  =  -[fo(sa)  -  r0(s)]  •  and  w2  =  — [rb(sfc)  -  f0(s)]  • 

ew.  For  a  source  position  s,  the  variables  uc  and  u'c  denote  the 
cone- beam  projections  of  xc  and  x'c  onto  the  detector  and  can 
be  obtained,  respectively,  by  replacing  x  with  xc  and  x'c  in 

=  +  l)  .  . 

w1(l-x)+w2(x  +  l)-  '  ’ 


The  rect  function  n c(x'c)  in  (8)  indicates  that  theMDFBP  al¬ 
gorithm  can  reconstruct  a  chord  image  from  knowledge  of  data 
only  on  the  cone- beam  projection  of  the  reconstruction  segment 
[xA, xb],  which  can  be  as  small  as  the  support  segment.  There¬ 
fore,  similar  to  the  BPF  algorithm,  the  M  DFBP  algorithm  can 
also  reconstruct  a  chord  image  from  data  containing  truncations 
[4],  [6], 


BPF  and  M  DFBP  algorithms  is:  1)  data  are  collected  over  the 
trajectory  segment  [sar»6]  and  2)  at  each  s,  data  only  on  the 
cone- beam  projection  of  the  reconstruction  segment  [xA,xB\  on 
the  chord  are  available.  It  follows  that,  because  the  reconstruc¬ 
tion  segment  [xA,xB\  can  be  chosen  as  small  as  the  support 
segment  [a;0l,xc2],  the  BPF  and  M  DFBP  algorithms  require, 
at  each  s,  data  only  on  the  cone-beam  projection  of  the  sup¬ 
port  segment  [xcl,xc2]  (instead  of  the  entire  chord-line  as  the 
chord-based  FBP  algorithm  requires).  Different selectionsof  the 
reconstruction  segment  [xA,xB]  imply  that  different  amounts  of 
data  at  each  s  can  be  used  for  reconstructing  the  chord  image. 
Under  the  ideal  continuous  conditions,  different  selections  of 
[. xA ,  xB\  yield  identical  chord  images.  However,  when  data  con¬ 
tain  noise  and  other  inconsistencies,  and  when  different  selec¬ 
tions  of  [xA,xB]  are  used,  the  BPF  and  M  DFBP  algorithms  in 
their  discrete  forms  may  yield  different  chord  images.  This  is  an 
issue  that  will  be  investigated  in  the  following. 

III.  Noise  Properties  of  Chord-Based 
Image  Reconstruction 

The  BPF,  M  DFBP,  and  FBP  algorithms  described  above  can 
be  applied  to  reconstructing  chord  images  from  parallel-,  fan-, 
and  cone-beam  data  [17],  A  Igorithms  analogous  to  the  B  PF  al¬ 
gorithm  that  are  capable  of  reconstructing  2- D  ROI  images  from 
truncation  data  have  also  previously  been  proposed  [8],  [13], 
[18].  We  study  the  noise  properties  of  chord-based  reconstruc¬ 
tion  by  use  of  these  algorithms  in  their  discrete  forms.  As  men¬ 
tioned  above,  the  BPF  and  M  DFBP  algorithms  can  reconstruct 
the  image  on  the  reconstruction  segment  [xA,xB]  as  long  as  it 
covers  the  support  segment  [xci,xc2].  We  analyze  image- noise 
properties  on  reconstruction  segments  of  different  lengths. 


C.  FBP  Algorithm 

The  chord-based  FBP  algorithm  [6],  [10]  can  be  expressed  as 


fc(Xci  Sa ,  $b) 


1 

2  7T2 


A 

I r-  f0(s)| 


uc  —  u'c  |rr  —  ro(s)|  dq 


D(r0(q)J) 


where  uc  indicates  the  cone-beam  projection  of  xc  onto  the  de¬ 
tector  and  is  determined  by  using  xc  to  replace  x  in  (8),  and  A 
denotes  the  di  stance  from  the  source  poi  nt  f0(s)  to  a  poi  nt  on  the 
detector  at  which  the  ray  connecting  r  and  f0(s)  intersects  the 
detector.  As  the  filtering  (i.e.,  the  integration  over  u'J  is  carried 
out  over  the  projection  of  the  straight  line  containing  the  chord, 
similar  to  other  existing  FB  P-based  algorithms,  thechord-based 
FBP  algorithm  cannot  exactly  reconstruct  ROI  images  from  data 
containing  transverse  truncations. 

D.  Data-Sufficiency  Conditions 

As  shown  in  (9),  a  data-sufficiency  condition  for  the  FBP 
algorithm  is:  1)  data  are  available  over  the  trajectory  segment 
s  e  [sa,  St]  and  2)  for  each  s,  data  on  the  cone- beam  projection 
of  the  chord  are  nontruncated.  This  condition  is  similar  to  that 
for  other  F BP-based  algorithms  [1],  [9],  [14],  [15],  [16],  From 
(4)  and  (7),  a  data-sufficiency  condition  for  the  chord-based 


A.  Analysis  of  Image-Noise  Properties 

Thechord-based  algorithms  invoke  three  major  mathematical 
operations:  differentiation,  backprojection,  and  filtration.  To  a 
large  extent,  the  BPF,  M  DFBP,  and  FBP  algorithms  differ  in  the 
orders  of  invoking  these  operations.  Below,  we  focus  on  inves¬ 
tigating  the  noise  properties  of  differentiation,  backprojection, 
and  filtration  in  the  BPF  algorithm.  The  approach  taken  in  the 
investigation  is  readily  applicable  to  analyzing  the  noise  prop¬ 
erties  of  the  M  DFBP  and  FBP  algorithms.  In  the  presence  of 
data  noise,  the  measured  projection  D(f0(s),/3)  should  be  in¬ 
terpreted  as  a  stochastic  process.  (Throughout  the  paper,  we  use 
boldfaceand  normal  letters  to  denote  a  stochastic  process  and  its 
mean,  respectively.)  Because  the  backprojection  g(xc,sa,sb ) 
and  the  final  image  fc(xc ,  sa,sb)  on  a  chord  are  computed  from 
D(r0(s),/3),  they  should  also  be  considered  as  stochastic  pro¬ 
cesses.  We  focus  in  this  section  on  investigating  the  chord-i  mage 
variance  fc(xc,  sa ,  sh)  by  the  investigation  of  noise  propagation 
through  each  step  involved  in  the  BPF  reconstruction  algorithm. 

1 )  Noise  Properties  of  Differentiation/Backprojection  for 
Parallel-Beam  Data:  Let  P (ud,s)  denote  the  parallel-beam 
projection  on  detector  bin  ud  acquired  at  view  s.  We  assume  the 
covariance  of  the  projection  data  P (ud,s)  to  be  uncorrelated, 
i.e., 

Cov  {P (ud,  s),  P  ( u'd ,  s')}  =  o2(ud,  s)S  ( ud  -  u'd )  S(s  -  s') 
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where  a2(ud,s)  denotes  the  variance  of  the  projection  data.  The 
backprojection  image  on  the  chord  is  given  by  [13] 


2 

f  d 

i,Sb)  =  J  ds— P(Md,  s) 


(ID 


where  ud  =  f  ■  eu(s).  The  final  image  variances  on  a  chord 
depend  upon  the  covariance  of  the  backprojection  image,  which, 
using  (11),  can  be  written  as 


Cov{g(a;c,  sa,  sb),g(x'c,  sa,  sj)}  =  0.  Therefore,  the 

covariance  of  the  backprojection  image  for  parallel-beam 
projections  can  be  re-expressed  as 

COV  {g(xc,  Sa,  Sb),  g  ( x'c,  Sa,  «&)}  ~  c(xc)8  (xc  -  x'c)  (16) 


where 


2 

c(xc)  =  aui  J  dsa2(ud,s). 


(17) 


COV  {g(xc,  Sa,  Sb),g  (x'c,  Sa,  Sfc)} 

7T  7T 

2  2 

=  J  ds  J  ds'CovlJ^P(ud,s),J^F(u'd,s')\.  (12) 

_  7T  7T  ^ 

2  2 

T he  eval uation  of  the  backprojection-i mage  covariance  i nvol ves 
the  data-derivative  covariance,  which  can  be  conveniently 
written  as 


2)  Noise  Properties  of  Differentiation/Backprojection  for 
Fan-Beam  Data:  In  the  fan-beam  case,  we  use  P (ud,s)  to 
denote  the  projection  on  detector  bin  ud  acquired  at  view  s. 
Again,  we  assume  P(«d,s)  to  be  uncorrelated  and  to  satisfy 
(10).  The  backprojection  image  in  (6)  can  be  re-expressed  as 

g(xc,sa,sb)  =  ds— — — jP '(ud,s)  (18) 

J  F-r0(s)| 


Cov 


AP(„JiS)i^pW,s0} 

xS  ( ud  -  u'd )  8(s  —  s')  +  T 


=  aui(j2(ud,  s ) 

para  id^d,  &  ) 


(13) 


where  ud  =  (Sr-  eM(s)/(f0(s)  -  r)  -  ew(s))  is  the  fan-beam 
projection  of  xc  onto  the  detector.  The  weighted-data  derivative 
p '(ud,s)  is  given  by 


where  Var{P(nd,  s)}  =  a2(ud,s)  denotes  the  known  data 
variance,  which  can  be  a  function  of  ud  and  s.  The  second 
term  Tpara(«d,  u'd,  s,  s')  represents  the  difference  between  the 
term  on  the  left-hand  side  and  the  first  term  on  the  right-hand 
side  of  (13).  A  Ithough  the  magnitude  of  Tpara(wd,  u'd,  s ,  s')  can 
be  larger  than  or  comparable  to  that  of  the  first  term  in  the 
right-hand  side  of  (13),  numerical  results  that  follow  show  that 
its  contribution  to  the  final  image  variance  on  a  chord  is  neg¬ 
ligibly  small.  Therefore,  we  consider  only  the  first  term  in  the 
derivation  of  the  chord-image  variance  below.  The  parameters 
a  and  cu  are  introduced  to  account  for  the  interpolation  effect 
of  the  discrete  data  derivative  and  discrete  backprojection  on 
the  chord-image  variance.  Using  the  first  term  in  (13),  we  can 
rewrite  (12)  as 


P'(ud,s)  =  A2(ud) 


df0(s) 

d 

P('UC^,  5) 

ds 

dud 

A(ud) 

(19) 


Using  (18),  one  can  write  the  covariance  of  the  backprojec¬ 
tion  image  as 


Cov  {g(xc,  Sa,Sb),  g  (x'c,  Sa,  S6)} 


Sb  Sb 


xCov{P  '(ud,s),P'(u'd,s')}  (20) 

whi ch  depends  upon  the  covariance  of  P'(ud,  s) .  A  gai n,  we  can 
conveniently  write  the  covariance  of  P '(ud,  s)  as 


COV  {g(xc,  Sa,  Sb),g  ( x'c ,  Sa,  Sfc)} 


«  au 


2 

J  dsa2(ud,  s)8  (ud  —  u'd) .  (14) 

_  7 r 
2 


Cov{P'(ud,s),P'(u'd,s')} 


=  acoa2(ud,s)A2(ud) 


df0(s) 


ds 


2 

6  (ud  -  u'd)  8(s  -  s') 


+  T/  an  (ud,  ud,  s,  s  ) 


(21) 


We  now  consider  two  points  xc  and  x'c  on  the  chord  and  let 
ud  and  u'd  denote  their  parallel-beam  or  fan-beam  projections 
onto  the  detector.  Clearly,  for  a  source  position  s  that  satisfies 
s  ^  sa  or  sb,  one  can  conclude  that 


ud  —  u'd  =  0  if  xc  —  x'c  =  0 

ud  u'd  ^  0  if  xc  -  x'c  ^  0.  (15) 


Thus,  if  xc  =  x'c,  Cov{g(a;c,  sa, 
8(0)aui  /  ds(j2(ud,  s)  and,  if 


s6),g(a;c,sa,sfc)} 

Xc 


where  Var{P(ttd,  s)}  =  a 2(ud,s)  denotes  the  known  data 
variance,  which  can  be  a  function  of  ud  and  s.  The  second 
term  T fan(ud,u'd,s,s')  represents  the  difference  between  the 
term  on  the  left-hand  side  and  the  first  term  in  the  right-hand 
side  of  (21).  As  numerical  results  indicate,  it  turns  out  that 
T fan(ud,u'd,s,s')  will  also  have  a  negligible  contribution  to 
the  chord-image  variance.  Therefore,  we  consider  only  the 
first  term  in  the  derivation  of  the  chord-image  variance  below. 
A  gai  n,  the  parameters  a  and  w  are  i  ntroduced  to  account  for  the 
^interpolation  effect  of  the  discrete  data  derivative  and  discrete 
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backprojection  on  the  chord-image  variance.  Using  the  first 
term  in  (21),  we  can  rewrite  (20)  as 


Cov  {g(xc,  Sa,  Sb),  g  (x'c,  Sa,  St)} 


f  ,  (T2(ud,s)A2(ud ) 
J  |r-f0(s)|4 


c(ro(s) 

ds 


2 

^  (iirf  —  lid)  • 


(22) 


Similar  to  the  parallel-beam  case  described  previously,  using 
(15),  one  can  conclude  that 


specific  interpolation  scheme  used  in  the  discrete  backprojec¬ 
tion.  We  illustrate  our  estimation  of  w  when  a  two-point  inter¬ 
polation  is  used  for  the  discrete  backprojection  in  the  parallel- 
and  fan-beam  cases.  One  can  readily  obtain  estimates  of  w  when 
other  interpolation  schemes  are  used.  In  our  studies,  ud  is  mea¬ 
sured  in  the  unit  of  detector-bin  size.  Let  *  =  floor(wd)  +  I0, 
where  J0  denotes  the  index  corresponding  to  ud  =  0.  At  a  back- 
projection  view  s,  we  use  P'  to  denote  the  discrete  weighed- 
data  derivatives.  For  a  given  ud,  we  express  the  interpolated 
weighted-data  derivative  as 


Cov  {g(xc,  sa,  Sb),  g  ( x'c ,  sa,  st)}  »  c(xc)6  (xc  -  x'c )  (23) 
where 


c(xc)  = 


AHud) 

I r-  f0(s)|4 


df0(s) 

ds 


2 

(J2(ud,s). 


(24) 


3)  Noise  Properties  of  Differentiation/Backprojection  for 
Cone-Beam  Data:  In  the  cone-beam  case,  let  P (ud,vd,s) 
denote  the  projection  at  view  s  on  a  detector  bin  specified  by 
(ud,vd).  In  the  so-called  “well  oriented”  rotation-coordinate 
system  [5],  unit  vector  eM(s)  is  parallel  to  and  unit  vectors 
ev(s)  and  ew(s)  are  orthogonal  to  the  tangential  direction 
(df0(s)/ds)  of  the  source  trajectory.  Let  ud  and  vd  denote 
the  coordinates  along  eM(s)  and  e„(s).  It  can  be  shown  [5] 
that  the  backprojection  image  depends  only  upon  the  data 
derivative  along  ud.  Therefore,  the  reconstruction  formula  for 
the  cone-beam  backprojection  image  can  be  obtained  from 
that  for  the  fan-beam  backprojection  image  in  (18)  by  simply 
replacing  P(ud ,  s)  and  A(ud)  with  P{ud, vd,  s)  and  A(ud,  vd), 
respectively.  Subsequently,  one  can  show  that  the  covariance  of 
the  cone-beam  backprojection  image  g(xc,sa,sb)  also  satisfies 
(23)  and  (24). 

4)  Estimation  of  Parameters  a  and  u  in  Discrete  Form:  The 
parameter  a  is  introduced  to  account  for  the  interpolation  ef¬ 
fect  of  the  discrete  data-derivative  on  the  chord-i  mage  variance. 
We  consider  a  two-point  derivative,  which  was  used  in  our  nu¬ 
merical  studies.  Let  P *  denote  the  discrete  data,  where  i  = 
1,2 and  I  indicates  the  total  number  of  data  points.  We 
assume  that  data  P,:  are  uncorrelated  and  use  Var{P;}  to  de¬ 
note  their  variances.  The  two-point  data  derivative  is  defined  as 

Pi^pi+i-Pi-i].  (25) 


Therefore,  the  variance  of  the  discrete  data  derivative  Pi  can  be 
written  as 


P'  =(1-7)P'+7P' 


+i 


(28) 


where  7  =  ud  -  %.  Furthermore,  we  can  write  the  variance  of 


as 


Var{P;s}  *  [(1  -  7)2  +  72]  Var{P'}.  (29) 

For  the  sake  of  simplifying  the  estimation  of  w,  we  have  ig¬ 
nored  the  correlation  between  P^  and  pf+1  and  assumed  that 
Var{P-}  «  Var{P'+1}.  We  selects  as  the  average  over  all  of 
the  possible  7s,  which  can  be  computed  as 

i 

w  =  J  [(1  -  7)2  +  72]  d7  =  (30) 

o 


Finally,  with  a  substitution  of  a  =  1/2  and  w  =  2/3  into  (17) 
and  (24),  we  obtain  the  variances  c(xc)  of  the  backprojection 
images  on  the  chords  for  the  parallel-beam  and  fan-beam  pro¬ 
jections,  respectively,  as 


2 

c)  =  \j  dsa2(ud,s) 


A^jud) 
r-  f0(s)|4 


dfpjs ) 
ds 


2 

a2(ud,s). 


(31) 

(32) 


5)  Noise  Property  of  Weighted  HilbertTransform  Over  a  Fi¬ 
nite  Interval:  The  weighted  H  i I bert transform  constitutes  an  im¬ 
portant  step  in  the  chord-based  BPF,  M  DFBP,  and  FBP  algo¬ 
rithms.  Consequently,  the  noise  properties  of  these  algorithms 
depend  upon  that  of  the  weighted  Hilbert  transform,  which  we 
study  in  the  following.  Let  f(xc,sa,sb)  denote  the  weighted 
H  i  I  bert  transform  of  the  backproj  ecti  on  i  mage  g(xc,  sa,sb) 


V„{P,}  =  lV»{P.+.}  +  Var{P,-.} 

When  data  variances  are  identical,  (26)  becomes 

Var{Pi}  =  ^Var{Pi}. 


(26) 


XB 

f{xc,sa,sb)  =  — 1 1—  (  -  b  Uv)  g  (a;/,  sa,  Sb ) .  (33) 

0[%c)  J  %c  Xc 

XA 


We  assume  that  g(x'c,sa,sb )  is  band-limited  to  vm.  Therefore, 
(27)  the  H  i  I  bert  transform  kernel  (l/a;c)  can  be  replaced  by 


Therefore,  in  our  studies,  we  select  a  =  1/2,  which  is  the  coef¬ 
ficient  in  (27). 

The  parameter  w  was  introduced  to  account  for  the  interpola¬ 
tion  effect  of  discrete  backprojection  on  the  chord-image  vari¬ 
ances.  The  estimated  value  of  u  depends  obviously  upon  the 


u m  '  2 

h(xc)  =  —irj  I  dusgn[u]e27r^Xc  =  (Tr^m^c)  ^  (34) 

J  xc 

~Pm 

In  the  presence  of  noise,  the  weighted  Hilbert  transform 
25f (xc,sa,sb)  should  be  interpreted  as  a  stochastic  process, 
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which  is  denoted  in  boldface.  The  variance  of  f (xc,sa,sh)  can 
be  written  as 


Var  5f,)|  =  b2^  y  j  dx'c  J  dr"b  (x'€)  b  (:r") 

X  h  (xc  —  x'c)  h  (xc  —  x” ) 

X  Cov{g(^,sa,s6),g(a;",sa,S6)}. 

(35) 

As  (16)  and  (23)  show,  the  backprojection  image  g(xc,sa,sb) 
can  be  treated  as  an  approximated  uncorrelated  stochastic 
process.  Using  the  result  of  (16)  or  (23),  we  can  write  (35)  as 


Var  |f(xc,  sa,s6)  j 


b2(xc) 


dx'cc  (x'c)  h2  (xc  —  x'c)  b 2  (a/c) .  (36) 


In  our  numerical  studies  in  Sections  1 1 1 -B 1  and  lll-Cl,  we  have 
used  vm  =  (i/2Ac),  where  Ac  denote  the  sample  interval  of 

S(xc,  sai  •%)• 

6)  Noise  Properties  of  Chord  Images:  Using  (4),  one  can 
write  the  variance  of  the  reconstructed  chord  image  as 


Var  {fc(xc,  s„,  St)} 
sa  Var  \f(xc,  sa,  sb  )]  + 


Var  {P(m^q,  vjg,  sa)} 
47 r2b2(xc) 


V  (/— xb)(xa~\~1)  !  v  (1-\-xa)(1-\-%b) 
x  JZ^c  +  •  lJ/J 


Substituting  (36)  into  (37)  yields 


Var  {fc(xc,  sa,  S5)} 

Xb 

J  dx'cc(x'c)  h2  (xc-x'c)  b2  (.<) 

VBl{P(UdO,VdO,Sa)} 

47 T2b2(xc) 

\/(1-xb)(xa+1)  ,  \/(1+xa)(1+xb)  nQ] 

— pt; — + — — J  (38) 

which  provides  a  formula  for  computing  the  chord-image 
variance. 


B.  Numerical  Studies  of  Noise  Properties  in  Parallel-Beam 
Reconstruction 

Using  the  parallel-beam  configuration  in  Fig.  1(b)  and  the  pa¬ 
rameters  given  in  Table  I,  we  calculated  noiseless  projections  for 
the  numerical  phantom  in  Fig.  1(a).  We  have  used  an  object-in- 
dependent  Gaussian  noise  model  and  an  object-dependent 
Poisson-noise  model  in  the  numerical  studies.  For  each  noise 
model,  we  generated  10  000  sets  of  noisy  data  by  using  noise¬ 
less  data  as  the  means.  The  standard  deviation  cr0  of  Gaussian 


Fig.  1.  (a)  Phantom  in  numerical  studies,  (b)  Parallel-beam  configuration,  (c) 
Fan-beam  configuration.  Solid  line  segment  with  endpoints  xA  and  xB  repre¬ 
sents  reconstruction  segment.  Thick  segment  with  endpoints  xcl  and  xc2  indi¬ 
cates  support  segment.  It  should  be  noted  that  [t  A .  <-„]  D  [a-d  Rectan¬ 
gular  ROI  is  decomposed  into  a  set  of  (dashed)  line  segments. 


whereas  the  standard  deviation  for  the  Poisson  noise  is  the 
noiseless  data  scaled  to  yield  a  total  photon  count  of  5  x  105  for 
each  view.  We  investigated  four  reconstruction  segments  with 
different  lengths  LAB  =  \xB  -  xA\  :  7.8, 10.0, 14.1,  and  20.0 
cm,  all  of  which  are  located  at  a;  =  4.06  cm.  It  can  be  observed 
in  Fig.  1(b)  that  the  length  of  the  support  segment,  5.5  cm  in 
length,  is  shorter  than  the  four  reconstruction  segments  consid¬ 
ered.  Therefore,  the  image  on  this  chord  can  be  reconstructed 
exactly  by  use  of  data  determined  by  these  reconstruction  seg¬ 
ments.  One  can  also  conclude  from  Fig.  1(b)  that  the  minimum 
data  required  by  the  first  three  reconstruction  segments,  which 
are  shorter  than  the  maximum  dimension  (about  15.6  cm)  of 
the  object  support,  are  truncated. 

1)  Noise  Properties  in  Reconstruction  from  Truncated  Par¬ 
allel-Beam  Data:  From  the  10  000  sets  of  data  containing 
Gaussian  noise,  we  used  (4)— (6)  to  reconstruct  10  000  noisy 
S5),  and  g(a^c,  5^,55),  respectively,  on 
the  four  reconstruction  segments  described  previously.  Based 
upon  these  noisy  reconstructions,  we  subsequently  computed 
their  corresponding  empirical  variances,  which  are  shown  in 
the  upper  row  of  Fig.  2.  We  compare  the  empirical  results  to 
the  analytical  results  obtained  by  use  of  (16),  (36),  and  (38). 
The  function  c(xc)  is  determined  by  using  a(ud,s)  =  a0  in 
(31),  where  a0  is  1.6%  of  the  maximum  value  in  noiseless 
data.  The  analytical  results  are  displayed  in  the  middle  row  of 
Fig.  2.  Similarly,  using  (4)— (6),  we  reconstructed  10  000  sets  of 
noisy  g(xc,  sa,  sb),  f sa,  sb),  and  fc(xm  sa,  sb)  on  the  four 
segments  from  10  000  sets  of  data  containing  Poisson  noise. 
The  computed  empirical  variances  from  these  noisy  images  are 
displayed  in  the  upper  row  of  Fig.  3.  Using  the  noiseless  data 
as  the  Poisson-noise  variance  a 2{ud,  s)  in  (17),  one  can  readily 
determine  c(xc);  and  using  the  determined  c(xc)  in  (16),  (36), 
and  (38),  one  can  compute  analytical  image  variances,  which 
are  displayed  in  the  middle  row  of  Fig.  3.  The  results  show 
that  the  analytical  and  empirical  results  agree  well  with  each 
other,  suggesting  that  (38)  provides  an  adequate  estimation  of 
the  chord-image  variance. 

It  can  also  be  observed  in  Fig.  2(c)  and  Fig.  3(c)  that  the 
shorter  the  reconstruction  segment,  the  higher  the  chord-image 
variances.  This  is  only  because  the  second  term  in  (38)  in¬ 
creases  as  Lab  (i.e.,  (xA-xc)(xc-xB))  decreases.  However, 
the  difference  of  the  chord-image  variances  in  the  central  part 
of  the  support  segment  is  quite  small  among  these  reconstruc¬ 
tion  segments.  The  implication  of  this  result  is  that  there  may 


noise  used  is  1.6%  of  the  maximum  value  in  the  noiseless  data,  26be  a  significant  gain  in  terms  of  dose  reduction  by  using  a  short 
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TABLE  I 

Parameters  for  Circular  Parallel-Beam,  Circular  Fan-Beam,  and  Helical  Cone-BEAM  Scanning  Configurations  Which  Are  Used 
in  our  Simulation  Studies.  SDD  Is  Source-to-Detector  Distance  and  SID  Denotes  Source-to-Isocenter  Distance 


pixel  position  (cm)  pixel  position  (cm)  pixel  position  (cm) 


a  b  c 

Fig.  2.  Empirical  (top  row)  and  analytical  (middle  row)  variances  of  (a)  g(a?c,  9a, sh),  (b)  f(xm &a,  sh),  and  (c)  fc(asc,3«,s&)  obtained  on  four  reconstruction 
segments  from  parallel-beam  data  containing  Gaussian  noise.  Difference  between  empirical  variances  and  analytical  variances  is  also  shown  in  bottom  row, 
which  demonstrates  analytical  variances  agree  well  with  the  empirical  variances.  Lengths  of  these  segments  are  indicated  in  the  box  in  upper-right  corners 
of  plots. 


reconstruction  segment,  because  data  required  to  reconstruct  an 
image  on  this  reconstruction  segment  is  less  than  that  required 
by  using  a  longer  reconstruction  segment,  thus  resulting  in  a 
reduced  illumination  coverage  to  the  object.  For  similar  X-ray 


intensities,  which  are  directly  related  to  the  data-noise  level,  the 
reconstruction  using  a  short  reconstruction  segment  appears  to 
yield  image  variance  within  the  support  segment  that  is  com¬ 
parable  to  that  obtained  with  a  longer  reconstruction  segment. 
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pixel  position  (cm)  pixel  position  (cm)  pixel  position  (cm) 


pixel  position  (cm) 


pixel  position  (cm) 


pixel  position  (cm) 


7.8  cm 
1CJ.0  cm 

14  T  cm 
20.0  cm 


0.1 


7  Bern 
10.0  cm 
-  ■  •  14.1  cm 
20.0  cm 


03 


7.8  cm 
10.0  cm 
14.1cm 
20.0  cm 


-1o 


pixel  position  (cm) 


ID 


-°4o 


pixel  position  (cm) 

b 


10 


-°4o 


pixel  position  (cm) 
c 


1C 


Fig.  3.  Empirical  (top  row)  and  analytical  (middle  row)  variances  of  (a)  g(zc,  ^}  ^}f  (b)  f(xc,  sa,  sb),  and  (c)  fc(xc,  sa,  sb)  obtained  on  four  reconstruction 
segments  from  parallel-beam  data  containing  Poisson  noise.  Difference  between  empirical  variances  and  analytical  variances  is  also  shown  in  bottom  row,  which 
demonstrates  analytical  variances  agree  well  with  empirical  variances.  Lengths  of  these  segments  are  indicated  in  the  box  in  upper-right  corners  of  plots. 


We  have  also  performed  numerical  studies  of  the  noise  prop¬ 
erties  of  the  reconstructed  R0I  images  by  use  of  the  BPF  and 
M  DFBP  algorithms  from  truncated  data.  Using  the  numerical 
phantom  in  Fig.  1  and  each  of  Gaussian-  and  Poisson-noise 
models  described  above,  we  generated  500  noisy,  truncated  data 
sets  for  image  reconstruction  on  reconstruction  segments  of  a 
length  LAb  =  10.0  cm,  as  shown  in  Fig.  1(a),  which  completely 
cover  the  ROI.  We  subsequently  reconstructed  500  noisy  images 
by  using  the  BPF  and  M  DFBP  algorithms.  Wedisplay  in  Fig.  4 
noisy  ROI  images  reconstructed  using  the  BPF  and  M  DFBP  al¬ 
gorithms  from  data  containing  Gaussian  noise  (upper  row)  and 
Poisson  noise  (lower  row). 

Using  the  reconstructed  500  sets  of  Gaussian-noise  images 
and  500  sets  of  Poisson-noise  images,  we  computed  empirical 
variances  of  the  ROI  images,  which  are  shown  in  the  upper  row 
and  lower  row  of  Fig.  5,  respectively.  We  display  in  the  third 
column  of  Fig.  5  the  variance  profiles  on  the  dashed  lines  indi¬ 
cated  in  the  variance  images.  Results  in  Fig.  5  support  the  con¬ 
clusion  that  both  BPF  and  M  DFBP  algorithms  yield  images  with 
comparable  variance  levels. 


a  b 


Fig.  4.  ROI -images obtained  by  useof  (a)  BPF  algorithm  and  (b)  M  DFBP  algo¬ 
rithm  from  truncated  parallel-beam  data  containing  Gaussian  noise  (upper  row) 
~8and  Poisson  noise  (lower  row). 
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Fig.  5.  Empirical  ROI-image  variances  obtained  by  use  of  (a)  BPF  and  (b)  M  DFBP  algorithms  from  truncated  parallel-beam  data  containing  Gaussian  noise 
(upper  row)  and  Poisson  noise  (lower  row),  respectively,  (c)  Variance  profiles  on  the  dashed  lines  indicated  in  columns  (a)  and  (b)  obtained  with  the  BPF  (solid) 
and  M  DFBP  (dashed)  algorithms,  respectively.  For  purpose  of  displaying  details  in  central  (i.e.,  low  variance)  regions,  we  have  applied  a  logarithmic  scale  to  the 
variance  images.  Display  windows  are  [-0.86,  0.50]  and  [-1.28,  0.20]  for  Gaussian  noise  and  Poisson  noise,  respectively. 


2)  Noise  Properties  in  Reconstruction  from  Nontruncated 
Parallel-Beam  Data:  As  discussed  above,  the  FBP  algorithm 
cannot  reconstruct  exactly  images  from  truncated  data.  There¬ 
fore,  we  study  the  noise  properties  of  the  FBP  algorithm  from 
parallel-beam  data  without  truncations.  For  the  purpose  of  com¬ 
parison,  we  have  also  included  reconstruction  results  of  the  BPF 
and  M  D  F  B  P  al  gorithms  from  the  same  nontruncated  data.  U  si  ng 
the  numerical  phantom  in  Fig.  1  and  each  of  the  Gaussian-  and 
Poisson-noise  models,  we  generated  500  noisy  data  sets  from 
which  500  noisy  images  were  obtained  by  use  of  each  of  the 
BPF,  M  DFBP,  and  FBP  algorithms.  Using  these  noisy  images, 
we  computed  empirical  variance  images,  which  are  shown  in  the 
upper  and  lower  rows  of  Fig.  6,  respectively,  for  the  Gaussian- 
and  Poisson-noise  models.  We  also  display  in  Fig.  7  the  variance 
profiles  on  the  dashed  lines  (i.e.,  on  a  chord)  indicated  on  the 
variance  images  in  Fig.  6.  The  profile  results  were  obtained  by 
use  of  the  BPF  (solid),  M  DFBP  (dashed),  and  FBP  (dotted)  al- 
gorithmsforthe(a)  Gaussian-noise model  and  (b)  Poisson-noise 
model .  1 1  can  be  observed  that  i  mage  vari  ances  obtai  ned  w  i  th  the 
three  algorithms  are  similar  and  that  the  only  difference  comes 
at  the  extreme  ends  of  the  shown  reconstruction  segments.  The 
BPF  and  M  DFBP  algorithms  show  a  significant  increase  in  the 
i mage  vari ance  at  both  ends  of  the  profile.  T he  reason  for  thi s  i s 
that  the  reconstruction  segment  was  taken  to  be  the  width  of  the 
image  array,  and  the  prefactor  for  the  finite  Hilbert  transform 
in  (4)  and  (7)  has  a  singularity  at  the  ends  of  the  reconstruction 
segment.  In  practical  situations  this  prefactor  is  of  little  con¬ 
sequence  because  the  reconstruction  segment  can  be  selected 
larger  to  avoid  the  singular  behavior;  furthermore,  because  the 
singularity  goes  as  the  -1/2  power,  its  effect  is  evident  only  very 
close  to  the  endpoints  of  the  reconstruction  segment. 


C.  Numerical  Studies  of  Noise  Properties  in  Fan-Beam 
Reconstruction 

Using  the  fan-beam  configuration  in  Fig.  1(c)  and  the  param¬ 
eters  listed  in  Table  I,  we  calculated  fan-beam  noiseless  data 
for  the  numerical  phantom  in  Fig.  1(a).  We  have  also  used  an 
obj ect-i  ndependent  G  aussi an-noi se  model  and  an  obj ect-depen- 
dent  Poisson-noise  model  in  this  numerical  study.  The  stan¬ 
dard  deviation  <r0  of  Gaussian  noise  used  is  2.3%  of  the  max¬ 
imum  value  in  noiseless  fan-beam  data,  whereas  the  standard 
deviation  for  the  Poisson  noise  is  the  noiseless  data  scaled  to 
yield  a  total  photon  count  of  5  x  105  for  each  view.  For  each 
noise  model,  10  000  sets  of  noisy  data  were  generated  by  use 
of  the  corresponding  noiseless  data  as  the  means.  We  investi¬ 
gated  reconstruction  segments  of  four  different  lengths  LAB  = 
7.8, 10.0, 14.1,  and20.0  cm.  A II  of  the  segments  are  located  at 
x  =  4.06  cm.  It  can  be  observed  in  Fig.  1(c)  that  the  length  of 
the  support  segment  is  5.5  cm,  which  is  shorter  than  the  four  re¬ 
construction  segments.  Therefore,  the  image  on  this  chord  can 
be  reconstructed  exactly  by  use  of  data  determined  by  these  re¬ 
construction  segments.  One  can  also  conclude  from  Fig.  1(c) 
that  data  determined  by  the  first  three  reconstruction  segments, 
which  are  shorter  than  the  maximum  di  mension  (about  15.6  cm) 
of  the  object  support,  are  truncated. 

1)  Noise  Properties  in  Reconstruction  from  Truncated  Fan- 
Beam  Data:  From  the  10000  sets  of  data  containing  Gaussian 
noise,  we  used  (4)— (6)  to  reconstruct  10  000  noisy  fc(xc,  sa,  sb), 
f(xc,  sa,  sb),  and  g(xc,  sa,  sb),  respectively,  on  the  four  recon¬ 
struction  segments  described.  Based  upon  these  noisy  recon¬ 
structions,  we  subsequently  computed  their  corresponding  em- 
29pirical  variances,  which  are  shown  in  the  upper  row  of  Fig.  8. 
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a  b  c 


Fig.  6.  Empirical  variance  images  obtained  by  use  of  (a)  BPF,  (b)  MDFBP,  and  (c)  FBP  algorithms  from  nontruncated  parallel-beam  data  containing  Gaussian 
noise  (upper  row)  and  Poisson  noise  (lower  row).  For  the  purpose  of  displaying  details  in  central  (i.e.,  low  variance)  regions,  we  have  applied  a  logarithmic  scale 
to  variance  images.  Display  windows  are  [-1.3,  0.30]  and  [-1.3,  0.40]  for  Gaussian  noise  and  Poisson  noise,  respectively. 


Fig.  7.  Variance  profiles  on  the  dashed  lines  indicated  in  variances  images 
shown  in  Fig.  6.  They  were  obtained  for  the  Gaussian  (a)  and  Poisson  (b)  noise 
models  by  use  of  the  BPF  (solid),  MDFBP  (dashed),  and  FBP  (dotted)  algo¬ 
rithms,  respectively. 


As  for  the  analytic  variance,  one  can  determine  c(xc )  by  using 
cr(ud,s )  =  (To  in  (32),  where  <r0  is  2.3%  of  the  maximum  value 
in  noiseless  fan-beam  data.  Using  c(xc )  in  (23),  (36),  and  (38), 
we  computed  analytically  image  variances,  which  are  displayed 
in  the  middle  row  of  Fig.  8.  Similarly,  using  (4)-(6),  we  recon¬ 
structed  10  000  sets  of  noisy  f c(xc,sa,sb),  f (xc,sa,sb),  and 
g(xc,  sa,  sb)  on  the  four  segments  from  10  000  sets  of  fan-beam 
data  containing  Poisson  noise.  The  computed  empirical  vari¬ 
ances  from  these  noisy  images  are  displayed  in  the  upper  row 
of  Fig.  9.  Furthermore,  using  the  noiseless  fan-beam  data  as  the 
Poisson-noise  variance  a2(ud,s)  in  (32),  one  can  readily  de¬ 
termine  c(xc).  Using  the  determined  c(xc )  in  (23),  (36),  and 
(38),  we  computed  analytically  image  variances,  which  are  dis¬ 
played  in  the  middle  row  of  Fig.  9.  It  can  be  observed  that 
the  analytic  and  empirical  results  agree  well  with  each  other, 


suggesting  that  (38)  provides  an  adequate  analytic  estimation 
of  the  chord-image  variance  for  the  fan-beam  case  as  well.  It 
is  interesting  to  note  in  Fig.  8(a)  and  Fig.  9(a)  that  the  vari¬ 
ances  of  g(xc ,  sa,  sb )  are  spatially  varying  on  the  chord.  Based 
upon  (32),  one  can  readily  conclude  that  this  spatial  variation  is 
caused  by  the  spatially  variant  factor  (A2(ud)/\f  -  r0(s)|4). 

Again,  from  these  results,  observations  similar  to  those  for 
the  parallel-beam  case  can  be  made  for  the  fan-beam  case.  For 
example,  as  Fig.  8(c)  and  Fig.  9(c)  show,  the  shorter  the  recon¬ 
struction  segment,  the  higher  the  chord-image  variances.  This 
is  only  because  the  second  term  in  (38)  increases  as  LAB  (i.e., 
(xA  -  xc)(xc  -  xB ))  decreases.  The  implication  of  these  re¬ 
sults  is  that  there  may  be  a  significant  gain  in  terms  of  dose  re¬ 
duction  by  using  a  short  reconstruction  segment,  because  data 
required  to  reconstruct  an  image  on  this  reconstruction  segment 
is  less  than  that  required  by  using  a  longer  reconstruction  seg¬ 
ment,  which  can  result  in  reduction  of  illumination  coverage  of 
the  object. 

We  have  also  performed  numerical  studies  of  the  noise  prop¬ 
erties  of  the  reconstructed  ROI  images  by  use  of  the  BPF  and 
MDFBP  algorithms  from  truncated  data.  Using  the  numerical 
phantom  in  Fig.  1  and  each  of  the  Gaussian-  and  Poisson-noise 
models  described  previously,  we  generated  500  noisy  truncated 
data  sets  for  image  reconstruction  on  reconstruction  segments 
of  a  length  LAB  =  10.0  cm,  as  shown  in  Fig.  1(a),  which  cover 
the  ROI  completely.  We  subsequently  reconstruct  500  noisy  im¬ 
ages  by  using  each  of  the  BPF  and  M  DFBP  algorithms.  Using 
the  reconstructed  500  sets  of  Gaussian-noise  images  and  500 
sets  of  Poisson-noise  images,  we  computed  empirical  variance 
images  within  the  ROI,  which  are  shown  in  the  upper  row  and 
lower  row  of  Fig.  10,  respectively.  M  oreover,  we  display  in  the 
30third  column  of  Fig.  10  the  variance  profiles  on  the  dashed  lines 
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Fig.  8.  Empirical  (top  row)  and  analytical  (middle  row)  variances  of  (a)  g(asc,sa,s&),  (b)  f(xc,  sa,  sb),  and  (c)  fc(xc,  saj  sb)  obtained  on  four  reconstruction 
segments  from  fan-beam  data  containing  Gaussian  noise.  Difference  between  empirical  variances  and  analytical  variances  is  also  shown  in  bottom  row,  which 
demonstrates  the  analytical  variances  agree  well  with  the  empirical  variances.  Lengths  of  these  segments  are  indicated  in  the  box  in  upper-right  corners  of  plots. 


indicated  in  the  variance  images.  Results  in  Fig.  10  support  the 
conclusion  that  both  BPF  and  M  DFBP  algorithms  yield  images 
with  comparable  variance  levels. 

2)  Noise  Properties  in  Reconstruction  From  Nontruncated 
Fan-BeamData:  The  FBP  algorithm  cannot  reconstruct  exactly 
images  from  truncated  fan-beam  data.  Thus,  we  evaluate  the 
noise  properties  of  the  FBP  reconstruction  from  nontruncated 
fan-beam  data.  For  the  purpose  of  comparison,  we  have  also 
included  reconstruction  results  of  the  BPF  and  M  DFBP  algo¬ 
rithms  from  the  same  data  sets.  Using  the  numerical  phantom 
and  fan-beam  configuration  described  above,  wegenerated  non¬ 
truncated  fan-beam  data  at  512  views  uniformly  covering  2w. 
Using  the  noiseless  data  as  the  means,  we  generated  500  sets  of 
data  containing  Gaussian  noise  and  500  sets  of  data  containing 
Poisson  noise.  The  standard  deviation  for  the  Gaussian  noise 
is  2.3%  of  the  maximum  value  of  the  noiseless  data,  whereas 
the  standard  deviation  for  the  Poisson  noise  is  the  noiseless 


data  scaled  to  yield  a  total  photon  count  of  5  x  105  for  each 
view.  For  a  given  chord  specified  by  sa  and  sb,  one  can  recon¬ 
struct  its  image  from  data  acquired  over  the  right-side  trajectory 
(i.e.,  s  g  [s0,s6]),  as  shown  in  Fig.  11(a).  Conversely,  one  can 
also  reconstruct  the  chord  image  from  data  acquired  with  both 
right-side  trajectory  (i.e.,  s  e  [sa,s6])  and  left-side  trajectory 
(i.e.,  s  g  [sfc,sa]),  as  shown  in  Fig.  11(b).  In  chord-based  image 
reconstruction,  we  decompose  image  area  into  chords  parallel 
to  the  vertical  direction  and  the  source  scans  from  sa  to  sb  and 
then  from  sb  to  sa,  as  shown  in  Fig.  11. 

For  each  chord  in  the  set  covering  the  image  area,  we  first  re¬ 
constructed  the  images  by  use  of  the  BPF,  M  DFBP,  and  FBP  al¬ 
gorithms  from  data  containing  Gaussian  noise  acquired  over  the 
right-side  trajectory  specified  by  s  g  [sa,s6].  Subsequently,  we 
computed  empirically  chord-image  variances  from  these  noisy 
reconstructions.  By  assembling  the  chord-image  variances,  we 
31  obtain  the  variance  images,  which  are  shown  in  the  upper  row  of 
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Fig.  9.  Empirical  (top  row)  and  analytical  (middle  row)  variances  of  (a)  g(a-c,  sa,  sb),  (b)  f(xc,  sa,sb),  and  (c)  tc(xc,  sa,sb)  obtained  on  four  reconstruction 
segments  from  fan-beam  data  containing  Poisson  noise.  Difference  between  empirical  variances  and  analytical  variances  is  also  shown  in  bottom  row,  which 
demonstrates  analytical  variances  agree  well  with  empirical  variances.  Lengths  of  these  segments  are  indicated  in  the  box  in  upper-right  corners  of  plots. 


Fig.  12,  for  the  B  PF,  M  DFBP  and  FBP  algorithms,  respectively. 
Similarly,  from  data  containing  Poisson  noise,  we  obtained  the 
image  variances,  which  are  displayed  in  lower  row  of  Fig.  12. 

We  show  in  column  one  of  Fig.  13  the  image  variances  on  a 
chord  specified  bysa  =  -7r/2and  sb  =  7r/2  obtained  from  data 
containing  Gaussian  noise  (upper  row)  and  Poisson  noise  (lower 
row),  respectively.  As  already  seen,  the  variance  increases  with 
the  position  along  the  chord  near  the  source  trajectory.  There 
is  little  difference  between  the  three  algorithms.  Furthermore, 
these  variance  images  have  similar  properties:  the  chords  on 
the  right  part  have  higher  and  more  nonuniform  image  vari¬ 
ance  than  those  on  the  left  part  in  the  image  area.  In  column 
2  of  Fig.  13,  we  show  the  profiles  on  the  middle  points  across 
the  vertical  chords  (i.e.,  on  the  middle  horizontal  lines  in  the 
variance  images  shown  in  Fig.  12).  The  results  reveal  that  some 
difference  of  the  M  DFBP  result  from  the  BPF  and  FBP  results 


in  the  peripheral  region.  This  difference  may  be  attributed  to 
the  data  weighting  prior  to  the  backprojection  step,  which  dif¬ 
fers  from  that  in  the  BPF  and  FBP  algorithms.  We  should  point 
out  that  this  difference  is  only  seen  in  the  extreme  periphery  of 
the  imaging  area.  For  most  practical  situations  these  three  algo¬ 
rithms  perform  virtually  identically  in  terms  of  image  variance. 

We  investigate  further  the  general  trend  of  the  variance  de¬ 
creasing  for  chords  on  the  left  of  the  variance  image.  The  first 
impression  of  this  behavior  is  that  this  trend  is  obvious,  be¬ 
cause  the  scanning  trajectory  is  on  the  right  side,  chords  on 
the  left  of  the  variance  image  are  reconstructed  with  a  longer 
scanning  trajectory.  It  appears  that  more  data  are  used  in  recon¬ 
structing  chords  covering  the  left  part  of  the  variance  image. 
This  explanation  is,  however,  incorrect.  First,  there  is  a  slight 
upturn  in  the  variance  for  the  chords  on  the  extreme  left  of 
32the  variance  image,  which  runs  counter  to  this  trend.  Second, 
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Fig.  10.  Empirical  ROI-image  variances  obtained  by  use  of  (a)  BPF  and  (b)  MDFBP  algorithms  from  truncated  fan-beam  data  containing  Gaussian  noise  (upper 
row)  and  Poisson  noise  (lower  row),  respectively,  (c)  Variance  profiles  on  dashed  lines  indicated  in  columns  (a)  and  (b)  obtained  with  BPF  (solid)  and  M  DFBP 
(dashed)  algorithms,  respectively.  For  the  purpose  of  displaying  details  in  central  (i.e.,  low  variance)  regions,  we  have  applied  a  logarithmic  scale  to  variance 
images.  Display  windows  are  [-0.76,  0.29]  and  [-0.75,  -0.04]  for  Gaussian  noise  and  Poisson  noise,  respectively. 


a  b 


Fig.  11.  (a)  Right-side  trajectory  (solid)  of  chord  (thick)  specified  bysa  and  sb. 
(b)  Right-side  trajectory  (solid)  and  I  eft- side  trajectory  (solid)  of  chord  (thick) 
specified  by  sa  and  sb.  Scanning  configuration  in  (b)  corresponds  to  a  full  fan- 
beam  scan. 


it  can  be  demonstrated  that  the  amount  of  data  going  into  the 
chord  reconstruction  does  not  necessarily  increase  as  the  scan¬ 
ning  trajectory  increases.  Based  upon  (32),  one  can  conclude 
that  the  true  cause  of  the  variance  behavior  is  spatially  depen¬ 
dent  weighting  factor,  l/|f -  f0(A)|,  in  the  BPF,  M  DFBP,  and 
FBP  algorithms  [19]. 

For  a  given  chord  specified  by  sa  and  sb,  when  full  scan  data 
are  available,  one  can  reconstruct  two  chord  images  by  use  of 
data  acquired  with  the  right-side  and  left-side  trajectories,  as 
shown  in  Fig.  11  and  then  obtain  a  final  chord  image  by  av¬ 
eraging  the  two  chord  images.  We  show,  in  Fig.  14,  the  variance 
images  of  the  full  scan  with  accompanying  profiles  in  Fig.  15. 

D.  Numerical  Studies  of  Noise  Properties  in  Cone-Beam 
Reconstruction 

The  BPF  and  MDFBP  algorithms  can  yield  exact  image 
reconstruction  on  a  chord  specified  by  sa  and  sb  as  long  as  the 


support  segment  on  the  chord  is  illuminated  by  the  X-ray  beam 
at  the  projection  views  s  e  [sa,si],  because  these  algorithms 
require  data  only  on  the  fan-beam  projections  of  the  support 
segment.  From  the  perspective  of  the  chord-based  algorithms, 
the  reconstruction  of  a  chord  image  from  cone-beam  data 
is  similar  to  that  of  a  chord  image  from  fan-beam  data.  In 
the  fan-beam  case,  the  orientation  of  the  fan-beam  planes  at 
different  views  remain  unchanged,  whereas,  in  the  cone-beam 
case  with  a  non-planar  trajectory,  the  orientation  of  the 
fan-beam-illumination  planes  generally  varies  from  view  to 
view.  As  discussed  in  Section  III-A3,  the  noise  properties  of 
differentiation,  backprojection,  and  filtration  in  the  cone-beam 
case  are  similar  to  that  in  the  fan-beam  case.  Therefore,  we 
include  below  only  the  study  results  on  the  noise  properties  of 
the  final  chord-images  reconstructed  from  cone-beam  data. 

1)  Helical  Cone-Beam  Configuration:  In  our  investiga¬ 
tion  of  the  noise  properties  of  image  reconstruction  from 
cone-beam  data,  we  consider  the  helical  trajectory,  which 
is  the  most  widely  used  in  clinical  and  industrial  CT.  For 
a  helical  scan,  the  source  trajectory  is  described  mathemat¬ 
ically  as  f0(s)  =  (i? cos s, i? sins, (h/2w)s),  where  R  is 
the  source  to  center-of-rotation  distance,  and  h  indicates 
the  helical  pitch.  For  a  chord  specified  by  sa  and  sb,  if 
(n  -  l)-7r  <  | -  sa|  <  (n  +  1 ) 7r ,  where  n  is  a  positive  odd 
integer,  the  chord  is  also  referred  to  as  an  mr-line  segment  [20], 
[21],  as  shown  in  Fig.  16(a).  In  particular,  when  n  =  1  and  thus 
0  <  |si-sa|  <  2-71- ,  the  chord  is  referred  to  as  a  %- line  segment 
[2],  [11],  In  this  paper,  we  consider  image  reconstruction  only 
on  7r-l i ne  segments  for  the  reason  that  the  imaging  volume 
enclosed  by  the  helix  can  be  filled  uniquely  and  completely  by 
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Fig.  12.  Empirical  variance  images  obtained  by  use  of  (a)  BPF,  (b)  M  DFBP,  and  (c)  FBP  algorithms  from  nontruncated  fan-beam  data  containing  Gaussian  (upper 
row)  and  Poisson  (lower  row)  noise.  For  the  purpose  of  displaying  details  in  central  (i.e.,  low  variance)  regions,  we  have  applied  a  logarithmic  scale  to  variance 
images.  Display  window  is  [-1.0,  0.65]  and  [-2.02,  0.29]  for  Gaussian  noise  and  Poisson  noise,  respectively. 
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Fig.  13.  Variance  profiles  along  central  (a)  vertical  and  (b)  horizontal  lines  in 
variance  images  shown  in  Fig.  12,  obtained  by  BPF  (solid),  M  DFBP  (dashed), 
and  FBP  (dotted)  algorithms  from  data  containing  Gaussian  noise  (upper  row) 
and  Poisson  noise  (lower  row). 


7r-line  segments  [11],  [12],  Thus,  7r-line  segments  can  be  used 
to  form  3-D  images  in  a  helical  cone-beam  scan. 

We  computed  the  noiseless  data  from  a  Shepp-Logan 
phantom  by  use  of  the  configuration  parameters  in  Table  I. 


Using  the  noiseless  data  as  the  means,  we  subsequently  gener¬ 
ated  500  sets  of  data  containing  Gaussian  noise  and  500  sets  of 
data  containing  Poisson  noise,  respectively.  The  standard  devi¬ 
ation  of  Gaussian  noise  is  chosen  to  be  0.7%  of  the  maximum 
value  of  the  noiseless  data,  whereas  the  standard  deviation  for 
the  Poisson  noise  is  the  noiseless  data  scaled  to  yield  a  total 
count  of  5  x  105  for  each  view. 

2)  Noise  Properties  in  Reconstruction  From  Helical  Cone- 
BeamData:  A  curved surfaceinthehelixvolumecanbeformed 
by  a  set  of  7r-line  segments  for  which  we  fix  one  end-point  at  sa 
and  sweep  the  other  endpoint  over  a  range  sb  e  [smi„,smax]. 
We  show  in  Fig.  16(a)  a  curved  surface  obtained  by  concate¬ 
nating  a  set  of  7r-line  segments  specified  by  sa  =  -i r  and 
sb  e  [— 0.5-71-, 0.57t].  Using  generated  noisy  helical  cone-beam 
data,  we  reconstructed  noisy  i mages  on  the  7r-line  surface  by  use 
of  the  BPF,  M  DFBP,  and  FBP  algorithms.  From  these  noisy  im¬ 
ages  we  subsequently  computed  empirical  image  variances  on 
the  7r-line  surface.  In  Fig.  17,  we  display  the  image  variances 
obtained  with  BPF,  M  DFBP,  and  FBP  algorithms  from  data  con¬ 
taining  Gaussian  noise  and  Poisson  noise. 

We  also  display  in  Fig.  18(a)  and  (b)  image  variances  on  the 
7r-line  segment  specified  by  sa  =  — k  and  sb  =  0  in  the  sur¬ 
face,  obtained  from  data  containing  Gaussian  noise  and  Poisson 
noi se,  respecti vel y .  T he  i  mage  vari ances  show  si  m i  I ar  character- 
istics  to  that  of  fan-beam  image  variances  observed  in  Fig.  13. 
Namely,  the  variance  image  on  the 7r-line  surface  in  Fig.  17  has 
a  structure  that  is  similar  to  the  right-side  scan  fan-beam  results 
presented  in  Section  lll-C-1;  the  images  on  7r-line  segments  re¬ 
constructed  from  smaller  helix  segments  tend  to  have  higher  and 
more  nonuniform  variances.  The  similarity  with  the  fan-beam 
34case  is  not  surprising  because  the  geometrical  arrangement  of 
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Fig.  14.  Variance  images  obtained  by  use  of  (a)  BPF,  (b)  M  DFBP,  and  (c)  FBP  algorithms  from  full-scan  fan-beam  data  containing  Gaussian  (upper  row)  and 
Poisson  (lower  row)  noise.  For  the  purpose  of  displaying  details  in  central  (i.e.,  low  variance)  regions,  we  have  applied  a  logarithmic  scale  to  variance  images. 
Display  windows  are  [-1.19,  0.31]  and  [-1.18,  0.20]  for  Gaussian  noise  and  Poisson  noise,  respectively. 
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Fig.  15.  Variance  profiles  along  the  central  (a)  vertical  and  (b)  horizontal 
lines  in  variance  images,  which  are  shown  in  Fig.  14,  obtained  by  BPF  (solid), 
M  DFBP  (dashed),  and  FBP  (dotted)  algorithms  from  data  containing  Gaussian 
noise  (upper  row)  and  Poisson  noise  (lower  row). 


the  7r-line  with  respect  to  its  scanning  trajectory  is  very  sim¬ 
ilar  to  the  relationship  between  the  chords  and  corresponding 
fan-beam  scanning  trajectory.  The  only  difference  is  that  there 
is  an  out-of-plane  bent  to  the  helix  segment. 


Fig.  16.  (a)  7r-line,  37r-line,  and  57r-linesegments  in  a  helical  scan,  (b)  Surfaces 
generated  in  imaging  volume  by  concatenating  jr-line  segments  specified  by 

sa  —  —  jr  and  s6  e  [ — 0 . 5 tt ,  0 . 5 tt] . 


Regarding  the  nonuniform  shape  of  the  variance,  one  can  at¬ 
tribute  the  high  variance  in  the  image  periphery  to  the  weighting 
factors  multiplying  the  data  derivatives  before  backprojection. 
As  the  algorithms  are  essentially  the  same  for  chord-image 
reconstruction  in  fan-  and  cone-beam  cases,  this  conclusion 
should  come  as  no  surprise.  In  the  2-D  fan-beam  case,  the 
variance  nonuniformity  and  level  was  reduced  by  equally 
weighting  reconstructions  for  both  left  and  right  side  scans  for 
each  chord  of  the  scanning  circle.  For  the  helical  configuration, 
it  is  clear  that  in  a  typical  scan  there  will  be  some  overscan  for 
nearly  all  the  chords  comprising  the  volume.  B ut  the  overscan 
part  of  the  trajectory  does  not  form  a  closed  loop  so  using  the 
overscan  data  to  reduce  image  variance  is  not  as  obvious  as 
the  case  of  the  circular  scan.  Future  work  will  focus  on  how 
to  utilize  the  overscan  data  for  non-closed  trajectories  for  the 
purpose  of  reducing  the  impact  of  data  noise  on  chord-base 
35ROI-image  reconstruction. 
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Fig.  17.  Empirical  variance  images  on  7r-line  surface  shown  in  Fig.  16(b)  obtained  by  use  of  (a)  BPF,  (b)  M  DFBP,  and  (c)  FBP  algorithms  from  data  containing 
Gaussian  noise  (upper  row)  and  Poisson  noise  (lower  row),  respectively.  For  the  purpose  of  displaying  details  in  central  (i.e.,  low  variance)  regions,  we  have  applied 
a  logarithmic  scale  to  variance  images.  Display  windows  are  [-1.80, 0.30]  and  [-1.80, 0.60]  for  Gaussian  noise  and  Poisson  noise,  respectively. 


Fig.  18.  Variances  on  central  vertical  line  segment,  specified  by  sa  =  -w 
and  sb  =  0,  obtained  with  BPF  (solid),  MDFBP  (dashed),  and  FBP  (dotted) 
algorithms  from  data  containing  (a)  Gaussian  noise  and  (b)  Poisson  noise. 


IV.  Discussion 

In  this  paper,  we  have  performed  analytic  and  numerical 
investigations  of  the  noise  properties  of  chord-based  image 
reconstructions  from  parallel-,  fan-,  and  cone-beam  data.  One 
of  the  main  points  of  the  investigation  was  to  test  whether 
or  not  the  reduced  illumination  in  designing  a  minimal  data 
set  for  a  particular  RO I  leads  to  a  significant  reduction  in 
exposure.  The  idea  was  to  compare  the  statistical  properties 
of  the  ROI  image  reconstructed  from  noise  realizations  of 
the  minimal  data  set  with  noise  realizations  of  the  full  data 
set.  Similar  noise  levels  were  used  in  both  data  sets,  which  are 
equivalent  to  modeling  similar  incidentX-ray  beam  intensities. 
Our  study  indicates  that  the  resulting  image  variance  was 
almost  the  same  for  images  reconstructed  from  both  data  sets. 


Thus,  the  minimal  data  set  for  ROI  reconstruction  leads  to  a 
significant  overall  dose  reduction,  because  the  body  is  exposed 
to  lower  amount  of  ionizing  radiation  in  the  reduced  scan. 
For  fan-beam  and  cone-beam  imaging,  we  explored  the  noise 
properties  of  the  extreme  periphery  of  the  imaging  region 
by  investigating  large  fan-  and  cone-angles.  Image  variance 
nonuniformity  was  found  to  be  caused  by  spatially  dependent 
weighting  factors  in  the  chord-based  reconstruction  algorithms. 
This  work  represents  a  study  of  the  noise  properties  of  chord- 
based  reconstruction  and  of  the  impact  of  physical  factors 
on  ROI  imaging  in  fan-beam  and  cone-beam  CT.  In  seeking 
ways  to  reduce  the  impact  of  noise  in  volume  imaging,  we  will 
investigate  schemes  to  incorporate  overscan  data.  The  analysis 
presented  in  this  work  can  directly  be  applied  to  chord-based 
image  reconstruction  for  general  trajectories.  Finally,  it  is 
important  to  investigate  the  behavior  of  the  ROI-reconstruction 
algorithms  when  other  important  factors  are  included  in  the 
data  model  such  as  X-ray  polychromaticity  and  nonlinear 
partial  volume  averaging. 
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Clinical  Relevance/ Application: 

Breast  tomosynthesis  has  received  renewed  interest  because  it  can  provide  3D 
information  about  the  breast.  This  work  concerns  iterative  reconstruction  of  accurate 
breast  tomosynthesis  images. 

Purpose: 

In  current  breast  tomosynthesis,  image  representation  with  non-isotropic  spatial 
resolution  is  used  for  reducing  computational  time.  This  can,  however,  lead  to  artifacts  in 
iterative  reconstruction  of  breast  tomosynthesis  images.  In  the  work,  we  investigate  the 
effect  of  non-isotropic  image  representation  on  the  reconstruction  accuracy.  Based  upon 
the  investigation,  we  devise  schemes  for  reducing  artifacts  in  iterative  reconstruction. 

Materials  and  Method: 

In  the  work,  we  focus  on  investigating  the  effect  of  non-isotropic  image  representation  on 
reconstruction  accuracy  of  iterative  algorithms.  The  iterative  algorithms  under  study 
include  the  total-variation  (TV)  based,  expectation  maximization  (EM),  and  algebraic 
reconstruction  technique  (ART)  algorithms.  Tomosynthesis  data  are  generated  at  12  and 
20  views  over  50  degrees  from  phantoms,  including  a  breast  phantom.  We  have 
reconstructed  images  by  using  image  representations  with  different  degrees  of  non¬ 
isotropic  spatial  resolution.  Specifically,  in  each  image  representation,  the  ratio  between 
the  in-plane  and  longitudinal  resolution  for  an  image  voxel  is  selected  to  be  a  value  less 
than  1. 

Results: 

We  have  reconstructed  images  by  use  of  TV-based,  EM,  and  ART  algorithms  for  image 
representations  with  different  ratios  of  in-plane  and  longitudinal  resolution.  Our  results 
demonstrate  that  non-isotropic  image  representation  can  lead  to  significant  artifacts  in 
reconstructed  images.  The  appearance  and  severity  of  the  artifacts  depend  not  only  upon 
the  ratio  between  the  in-plane  and  longitudinal  resolution  but  also  upon  the  iterative 
algorithms.  The  TV-based  algorithm  seems  to  be  less  susceptible  to  the  effect  than  the 
EM  and  ART  algorithms.  Through  the  selection  of  algorithm  parameters,  the  artifacts  can 
be  reduced. 

Conclusion: 

The  non-isotropic  image  representation  can  significantly  affect  reconstruction  accuracy 
obtained  with  iterative  algorithms  in  breast  tomosynthesis. 
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Abstract — Recently,  a  reduced  circular  sinusoid  scan  was  proposed 
in  which  data  are  collected  only  over  an  open  portion  of  the 
full  circular  sinusoid  trajectory.  In  general,  for  reconstructing  the 
same  ROI  image,  the  reduced-scan  approach  uses  less  cone-beam 
data  than  what  required  for  a  full  circular  sinusoid  scan.  It  is  of 
practical  significance  to  investigate  and  evaluate  the  numerical  and 
noise  properties  of  the  image  reconstruction  for  reduced  circular 
sinusoid  trajectories.  In  this  work,  we  carry  out  studies  on  the 
noise  properties  in  images  for  full  and  reduced  circular  sinusoid 
trajectories  reconstructed  by  use  of  the  recently  developed  full- 
scan  and  reduced-scan  FBP  algorithms.  It  can  be  observed  that 
the  spatially  dependent  weighting  factors  in  the  FBP  algorithms  can 
lead  to  non-uniform  image  variances. 


I.  Introduction 

Circular  trajectory  cannot  yield  data  sufficient  for  exact  re¬ 
construction  of  3D  images.  On  the  other  hand,  non-circular 
trajectories  such  as  helical  and  saddle  trajectories  can  produce 
data  that  would  allow  for  exact  3D  image  reconstruction.  In  the 
last  several  years,  there  have  been  significant  advances  in  the 
development  of  algorithms  for  image  reconstruction  from  data 
acquired  with  non-circular  trajectories  [1],  [2],  [3].  It  is  expected 
that,  among  these  non-circular  trajectories,  the  class  of  circular 
sinusoidal  trajectories,  which  includes  the  saddle  trajectory  as  a 
special  case,  may  find  important  applications  in  cardiac  imaging 
and  C-arm  CT.  Based  upon  the  Pack-Noo  formula  [4],  [5],  Yang 
et.  al.  have  developed  a  filtered  backprojection  (FBP)  algorithm 
for  reconstructing  images  within  regions  enclosed  by  the  circular 
sinusoid  trajectory,  which  we  refer  to  as  a  full  scan  [6].  In  the  last 
year,  we  proposed  a  reduced  circular  sinusoid  scan  in  which  data 
are  collected  only  over  an  open  portion  of  the  circular  sinusoid 
trajectory.  Using  the  result  in  [5],  [6],  we  have  subsequently  de¬ 
rived  a  reconstruction  formula  for  achieving  exact  reconstruction 
of  images  in  a  reduced  scan.  In  general,  for  reconstructing  the 
same  ROI  image,  our  reduced-scan  approach  uses  less  cone-beam 
data  than  what  required  for  a  full  circular  sinusoid  scan  [6], 
thus  leading  a  reduction  of  radiation  dose  to  the  imaged  subject, 
scanning  time,  and  possibly  imaging  artifacts.  It  is  of  theoretical 
and  practical  significance  to  investigate  and  evaluate  the  numerical 
and  noise  properties  of  the  algorithms  for  full  and  reduced  circular 
sinusoid  trajectories.  In  the  absence  of  data  inconsistencies  such 
as  noise  and  sample  aliasing,  the  algorithms  for  data  acquired 
with  full  and  reduced  circular  sinusoid  scans  yield  identical 
reconstruction.  However,  in  the  presence  of  data  inconsistencies, 
the  algorithms  for  full  and  reduced  circular  sinusoid  trajectories 
are  expected  to  behave  differently,  thus  yielding  images  with 
different  noise  properties.  In  this  work,  we  have  carried  out 
studies  on  the  noise  properties  in  images  for  full  and  reduced 
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Figure  1.  Circular  sinusoidal  trajectories  of  order  m  =  3.  (a)  A  full-scan  (or, 
equivalently,  closed)  circular  sinusoidal  trajectory,  where  A  6  [— 7r,7r].  (b)  A 
reduced-scan  (or,  equivalently,  open)  circular  sinusoidal  trajectory,  where  A  E 
[Aa,  At,],  and  — 7r  <  Aa  <  A b  <  n. 


circular  sinusoid  trajectories  reconstructed  by  use  of  the  recently 
developed  full-scan  and  reduced-scan  FBP  algorithms. 

II.  Circular  sinusoidal  trajectory 

In  a  fixed-coordinate  system,  a  circular  sinusoidal  source  tra¬ 
jectory  ro(A)  of  order  m  can  be  expressed  mathematically  as 

ro(A)  =  (.RcosA,  .RsinA,  focosmA)T,  (1) 

where  A  denotes  the  rotation  angle  of  the  source  point,  R  indicates 
the  distance  from  the  source  to  the  rotation  axis,  h  is  the  amplitude 
of  the  source’s  oscillation  along  2  direction,  and  the  superscript 
T  stands  for  a  matrix  transpose.  A  circular  sinusoidal  trajectory 
has  a  minimal  period  2i r,  and,  as  shown  in  Fig.  la,  a  closed 
circular  sinusoidal  trajectory  is  formed  when  the  rotation  angle  A 
varies  over  the  entire  period  of  2i r,  which  we  refer  to  as  a  full- 
scan  circular  sinusoidal  trajectory.  In  contrast,  an  open  circular 
sinusoidal  trajectory  is  formed,  as  shown  in  Fig.  lb,  when  the 
scanning  angular  range  is  smaller  than  2i r,  which  we  refer  to  as 
a  reduced- scan  circular  sinusoidal  trajectory.  Notice  that,  when 
m  =  2,  the  circular  sinusoidal  trajectory  becomes  the  standard 
saddle  trajectory  [6],  [7]  and  that,  when  h  becomes  zero,  a  circular 
sinusoidal  trajectory  becomes  the  2D  circular  trajectory. 

The  cone-beam  projection  of  an  object  function  f(f)  at  rotation 
angle  A  is  defined  as 

POO 

g(X,0)  =  dtf(ro(X)+t0),  (2) 

Jo 

where  6  E  S2,  and  S 2  is  a  set  of  all  unit  vectors  over  a  hemisphere 
in  3D  space. 
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III.  The  reconstruction  algorithm  for  a  full  scan 

AND  A  REDUCED  SCAN 


Consider  a  full  scan  with  a  closed  circular  sinusoidal  trajectory 
as  shown  in  the  left  panel  of  Fig.  2.  A  plane  II  intersects  with 
the  trajectory  at  n(>  3)  different  points  ro(A^),  i  =  1 ,  ...,n,  and 
Ai  <  A2...  <  An.  A  convex  polygon  Un  is  formed  by  these 
intersection  points  within  the  plane.  It  has  been  shown  [6]  that 
image  on  this  polygon  can  be  reconstructed  as 

1  1 

f(x)  =  A  }C(x,  e»,  Ai;  Ai+i)  +  -K{x,  en,  A„,  2n  +  Ai),  (3) 
z  i= 1 

where 


K{x,e,\~  ,\+)  ^  ~ 


A-  CZA||^-r0(A)|| 


gF{ \,x,e),  (4) 


and  the  gp(\,x,e)  denotes  the  filtered  projection  data,  given  by 

1  dg{\,0(\,x,ei,'))) 


gF(\,x,ii)  =  J 


d~f- 


smy 


d\ 


(5) 


0(X,x,  (7,7)  =  COS7  a(X,  x)  +  sin7/3(A,  x,  e*), 


and 


/5(A,  x,  e$) 


&(A,  a?) 


—  (e$  •  a(X,x))a(X,x) 
\e$  -  (ei  •  a(A,z))d(A,z)| 


r0(A) 


||^-ro(A)||' 


(6) 

(7) 

(8) 
i  = 


The  filtering  directions  are  given  by  e*  =  ||£°(a-+|)7o(a*)|| 

1  77  —  1  and  e  — 

i, n  i  ana  en  -  ||f0(27r+A1)-ro(A„)||  • 

However,  for  a  reduced  scan,  data  cannot  be  collected  over 
the  trajectory  segment  between  An  and  2tt  +  Ai,  as  shown  in 
the  right  panel  of  Fig.  2.  Although  the  actual  trajectory  segment 
between  An  and  27 r  +  Ai  is  missing,  it  can  be  shown  that  the 
contribution  of  the  data  on  the  trajectory  segment  over  the  angular 
range  [An,  2tt  +  Ai]  can  be  replaced  by  that  of  the  open  trajectory 
segment  over  angular  range  [Ai,  An]  (i.e.  JC(x,  en,  An,  2tt-\-  Ai)  = 
K{x,  en,  Ai,  An)).  Based  upon  this  observation,  we  can  obtain  the 
reduced-scan  FBP  reconstruction  algorithm  as 


n— 1 


/ (**0  —  r)  ^  ^  JC(x,  ei,  A 2,  Ai+i)  -|-  K,(x,  en ,  Ai ,  An) 


(9) 


i=l 


IV.  Numerical  results 

We  have  performed  numerical  studies  to  investigate  the  noise 
properties  of  the  full-scan  and  reduced-scan  FBP  algorithms.  In 
these  studies,  a  low-contrast  Shepp-Logan  phantom  was  used  to 
generate  the  noiseless  cone-beam  projection  data  for  the  circular 
sinusoidal  trajectories  of  orders  m=2  (i.e.,  saddle  trajectory)  and 
m=3.  We  generated  the  full-scan  data  at  300  projection  views 
uniformly  distributed  over  [—7 r,  ir\.  A  portion  of  the  full-scan  data 
over  the  angular  range  [—  jtt,  jtt]  was  used  to  simulate  the  data 
acquired  in  a  reduced  scan.  By  adding  stationary  Gaussian  noise 
to  the  noiseless  data,  we  produced  500  sets  of  noisy  data  for 
each  scan.  From  these  full- scan  and  reduced- scan  data  sets,  we 


Figure  2.  Illustration  of  the  source  trajectories  in  a  full  scan  (left)  and  a  reduced 
scan  (right)  for  a  circular  sinusoidal  trajectory  of  orders  m=2  (top)  and  m=3 
(bottom).  The  plane  £2nis  parallel  to  the  x  —  y  plane,  and  intersects  with  the 
source  trajectory  at  four  (m=2)  or  six  (m=3)  different  points,  ro(A?;),  i  =  [1  ...n]. 
The  vectors  e*,  %  =  [l...n],  denote  the  filtering  direction. 


reconstructed  images  by  use  of  the  FBP  algorithm  described  in 
Sec.  III. 

In  Fig.  3,  we  displayed  one  of  the  reconstructed  images  within 
a  slice  at  2  =  0  cm  for  the  full- scan  and  reduced- scan  circular 
sinusoidal  trajectories  of  order  m=2.  Using  the  images  recon¬ 
structed  from  these  noisy  data  sets,  we  have  computed  empirical 
image  variances,  which  are  shown  in  Fig.  4.  For  the  quantitative 
comparison,  the  image  variances  along  the  middle  horizontal  line 
and  the  middle  vertical  line  are  also  displayed  in  Fig.  4.  It  can  be 
observed  in  Fig.  4,  that  the  image  variances  obtained  from  the  full- 
scan  data  and  reduced-scan  data  by  use  of  the  FBP  algorithms  are 
non-uniform.  The  reason  for  this  phenomena  is  that  the  spatially 
dependent  weighting  factors  (i.e.,  l/\x  —  fo(X)\)  involved  in  the 
FBP  algorithms  can  amplify  data  noise  and  inconsistencies  [8], 
[9].  One  can  also  observe  that,  for  the  full  scan,  the  empirical 
variance  within  the  central  region  of  the  field  of  the  view  (FOV) 
is  relative  smaller  compared  to  that  within  the  peripheral  region, 
and  the  variance  increases  as  the  distance  to  the  center  of  rotation 
increases.  However,  for  a  reduced  scan  in  which  the  trajectory 
is  on  the  right-hand  side  of  the  FOV,  the  spatially  dependent 
weighting  factors  results  in  image  variance  on  the  left  side  is 
lower  than  that  on  the  right  side  of  the  image  space. 

Again,  we  showed  in  Fig.  5  the  reconstructed  images  within 
a  slice  at  2  =  0  cm  for  the  circular  sinusoidal  trajectory  of 
order  m=3.  The  empirical  image  variances  and  the  corresponding 
profiles  are  displayed  in  Fig.  6.  The  similar  behavior  of  the  image 
variances  can  be  observed. 

V.  Discussion 

In  this  work,  we  have  investigated  the  noise  properties  in 
images  reconstructed  by  use  of  the  FBP  algorithms  for  full  and 
reduced  circular  sinusoid  scans.  It  was  observed  that  the  spatially 
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Figure  3.  The  reconstructed  images  within  the  2D  slices  at  z  =  0  cm  from 
the  noisy  projection  data  acquired  in  a  full-scan  (left)  and  a  reduced-scan  circular 
sinusoidal  trajectory  of  order  m=2  (right).  The  display  grey  scale  is  [1.0,  1.04]. 


Figure  6.  Top  row:  image  variances  within  the  2D  slices  at  z  —  0  cm  obtain 
from  the  noisy  projection  data  in  a  full  scan  (left)  and  a  reduced  scan  (right). 
Bottom  row:  Bottom  row:  variance  profiles  on  the  middle  horizontal  line  (left) 
and  the  middle  vertical  line  (right)  of  the  variance  images  in  a  full  scan  (dashed 
line)  and  a  reduced  scan  (solid  line)  for  the  circular  sinusoidal  trajectory  with 
order  m=3. 


dependent  weighting  factors  in  the  FBP  algorithms  can  lead  to 
non-uniform  image  variances. 
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Figure  4.  Top  row:  image  variances  within  the  2D  slices  at  z  —  0  cm  obtain 
from  the  noisy  projection  data  in  a  full  scan  (left)  and  a  reduced  scan  (right). 
Bottom  row:  Bottom  row:  variance  profiles  on  the  middle  horizontal  line  (left) 
and  the  middle  vertical  line  (right)  of  the  variance  images  in  a  full  scan  (dashed 
line)  and  a  reduced  scan  (solid  line)  for  the  circular  sinusoidal  trajectory  with 
order  m=2. 


Figure  5.  The  reconstructed  images  within  the  2D  slices  at  z  =  0  cm  from 
the  noisy  projection  data  acquired  in  a  full-scan  (left)  and  a  reduced-scan  circular 
sinusoidal  trajectory  of  order  m=3  (right).  The  display  grey  scale  is  [1.0,  1.04]. 
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Abstract 

In  classical  tomosynthesis,  the  X-ray  source  generally  is  moved  along  a  curve  segment,  such  as  a 
circular  trajectory,  within  a  plane  that  is  perpendicular  to  the  detector  plane.  Studies  suggest  that, 
when  the  angular  coverage  and  number  of  projection  views  are  limited,  it  can  be  difficult  to  recon¬ 
struct  accurate  images  within  planes  perpendicular  to  the  detector  plane  in  classical  tomosynthesis.  In 
this  work,  we  investigate  imaging  strategies  in  tomosynthesis  using  trajectories  that  are  not  confined 
within  a  plane  perpendicular  to  the  detector  plane.  We  expect  that  such  trajectories  can  increase 
data  information  and  thus  lead  reconstructed  images  with  improved  quality.  Numerical  studies  were 
conducted  for  evaluating  the  image-reconstruction  quality  in  classical  tomosynthesis  and  tomosynthesis 
with  trajectories  that  are  not  confined  within  a  plane  perpendicular  to  the  detector  plane.  The  results 
of  the  studies  indicate  that,  with  the  same  number  of  views,  (or,  equivalently,  the  same  amount  of 
imaging  radiation),  data  acquired  in  tomosynthesis  with  the  trajectories  that  are  not  confined  within  a 
plane  perpendicular  to  the  detector  plane  generally  contain  more  information  than  that  acquired  with 
classical  tomosynthesis  and  can  thus  yield  images  with  improved  quality. 


1  Introduction 

In  the  past  decade  or  so,  there  have  been  renewed  interests  in  the  development  of  tomosynthesis  for  appli¬ 
cations  to  breast  imaging,  image-guided  radiation  therapy,  and  security  scans  [1,  2].  In  these  studies,  the 
X-ray  source  often  is  moved  along  a  curve  segment  within  a  plane  that  is  perpendicular  to  the  detector 
plane.  For  convenience,  we  refer  to  such  a  source  trajectory  as  a  planar  trajectory.  In  some  practical 
applications  of  tomosynthesis,  however,  the  source  trajectory  need  not  be  restricted  to  be  within  a  planar 
trajectory,  and  it  can  be  a  trajectory  in  3D  space.  For  example,  the  X-ray  source  may  move  over  a  sur¬ 
face  such  as  a  portion  of  the  sphere  in  some  applications.  In  these  cases,  we  refer  to  the  trajectories  as 
non-planar  trajectories.  We  hypothesize  that  tomosynthesis  with  certain  non-planar  trajectories  can  yield 
images  with  improved  quality  over  that  obtained  with  a  planar  trajectory.  In  this  work,  we  conduct  a  pre¬ 
liminary  investigation  of  image  reconstruction  in  classical  tomosynthesis  that  uses  a  planar  trajectory  and 
tomosynthesis  that  uses  a  non-planar  trajectory.  We  will  apply  the  TV-minimization  iterative  algorithm 
for  image  reconstruction  from  data  acquired  in  the  two  tomosynthesis  configurations  [3,  4]. 

2  Scanning  configurations  for  tomosynthesis 

In  classical  tomosynthesis,  the  x-ray  source  is  moved  along  a  curve  segment.  When  the  curve  segment  is  a 
portion  of  a  circular  trajectory,  it  can  be  expressed  mathematically  as 

ro(A)  =  (R sinA,  0,i?cosA)  A  G  [As,Ae],  (1) 

where  A  denotes  the  rotation  angle,  R  is  the  distance  from  the  source  point  to  the  rotation  center,  and 
parameters  As  and  Ae  denote  the  starting  and  ending  angles  of  the  rotation.  For  each  rotation  angle  A, 
a  2D  detector  is  used  to  measure  the  projection  data,  which  is  the  path  integral  of  the  x-ray  attenuation 
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coefficient  along  the  rays  connecting  the  source  spot  and  an  individual  detector  bin.  In  terms  of  discrete 
detector  array  and  image  array,  the  data  model  can  be  written  as  the  weighted  sum  over  the  pixels  traversed 
by  the  source-bin  ray,  i.e., 

9  =  Mf,  (2) 

where  /  and  g  are  two  finite  vectors,  representing  the  image  function  and  projection  data.  The  dimensions 
of  /  and  g  are  7Vimage  and  iVdata,  which  denote  the  elements  of  discrete  image  and  data  arrays.  The  system 
matrix  M  thus  composes  iVdata  row  vectors  Mi  of  dimension  of  7Vimage,  yielding  gi  =  Mi  •  /.  Each  element 
Mij  of  the  system  matrix  is  the  length  of  the  ith  ray  traversing  the  jth  pixel.  The  process  of  image 
reconstruction  is  to  obtain  an  image  represented  by  the  finite  vector  /  from  knowledge  of  the  data  vector 
g  and  the  system  matrix  M . 

In  classical  tomosynthesis  for  breast  imaging,  projection  data  are  measured  at  about  20  projection  views 
uniformly  distributed  over  an  angular  range  of  about  30°.  One  may  improve  the  image  quality  of  classical 
tomosynthesis  by  increasing  the  number  of  projections  and/or  the  scanning  angular  range.  However,  this 
approach  results  in  increased  imaging  radiation  dose  to  the  subject.  On  the  other  hand,  it  may  be  possible 
to  design  innovative  scanning  configurations  for  acquiring  more  information  than  that  collected  with  a 
planar  trajectory  in  Fig.  la  in  classical  tomosynthesis,  yet  without  increasing  the  imaging  radiation  dose. 
In  this  work,  we  investigate  a  non-planar  imaging  configurations  in  tomosynthesis  for  acquiring  increased 
data  information  and  thus  for  achieving  improved  image  quality.  In  contrast  to  the  planar  trajectory  in 
classical  tomosynthesis,  we  study  in  this  work  a  non-planar  trajectory  that  is  on  a  curved  surface.  In 
particular,  we  consider  a  source  trajectory  that  consists  of  two  orthogonal  curve  segments  on  a  portion  of  a 
spherical  surface,  as  shown  in  Fig.  lb.  The  two  curve  segments  of  the  non-planar  trajectory  are  expressed 
as 

ro(A)  =  (.RsinA,  0,  R cosA)  A  G  [As,Ai],  (3) 

r0( A)  =  (0,  i?sin(A),  Rcos(X))  A  G  [A2,Ae].  (4) 

It  can  be  observed  that  the  two  circular  segments  are  chosen  to  be  within  x-z  and  y-z  planes.  Both  of 
them  are  on  a  spherical  surface  of  radius  R.  Clearly,  they  are  not  within  a  single  plane.  In  this  work, 
we  will  investigate  and  compare  images  reconstructed  by  use  of  the  TV-minimization  algorithm  from  data 
acquired  with  the  planar  trajectory  in  classical  tomosynthesis  and  from  data  obtained  with  the  non-planar 
trajectory. 

3  Reconstruction  algorithm 

We  briefly  review  the  TV-minimization  algorithm  [3]  that  is  used  to  reconstruct  images  from  tomosynthesis 
data  in  this  work.  In  many  imaging  applications,  it  is  not  uncommon  that  images  to  be  reconstructed  are 
relatively  constant  over  extended  volumes,  and  that  significant,  rapid  variation  in  the  image  may  occur 
only  at  boundaries  of  internal  structures.  In  these  situations,  the  image  formed  by  taking  the  magnitude  of 
its  gradient  could  be  approximately  sparse.  Therefore,  the  reconstruction  strategy  considered  in  the  TV- 
minimization  algorithm  is  to  incorporate  the  assumption  of  gradient  image  sparseness  on  the  image  function 
/  to  arrive  at  a  solution  from  knowledge  of  the  data  g.  Based  upon  this  strategy,  the  TV-minimization 
algorithm  seeks  to  find  the  solution  for  the  optimization  problem  below: 

/*  =  argmin||/||Ty,  (5) 


with  two  constraints, 

\Mf-g\<e  and  /»  >  0, 

where  /*  is  the  reconstructed  image.  The  inequality  used  in  the  first  constraint  accounts  for  data  inconsis¬ 
tency,  such  as  noise,  continuous-to-discrete  inconsistency.  The  parameter  e  can  be  selected  for  controlling 
the  impact  level  of  potential  data  inconsistency  on  the  image  reconstruction.  In  our  investigation  of  3D 
tomosynthesis  reconstruction,  the  TV  of  a  discrete  image  is  defined  as 


WfWrv 


E  |V/a,t,r|  =  E 


fs— l,t,r)2  +  ( fs,t,r  fs,t— l,r)2  +  ( fs,t,r  fs,t,r— l)2? 


(6) 
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Figure  1:  (a)  Scanning  configuration  with  a  planar  trajectory  within  x-z  plane  in  classical  tomosynthesis, 
(b)  Scanning  configuration  with  a  non-planar  trajectory  that  consists  of  two  orthogonal,  circular  segments 
within  x-z  and  y-z  planes,  respectively  (right). 


where  8,  £,  and  r  indicate  the  pixel  indices  within  the  3D  image. 

The  implementation  of  the  TV-minimization  algorithm  includes  two  major  steps:  gradient  descent 
method  and  projection  on  convex  sets  (POCS)  (see  e.g.  Sec.  15.4.5  of  Ref.  [5]).  The  gradient  descent 
method  is  used  for  minimizing  the  image  total  variation,  whereas  the  POCS  is  used  for  enforcing  the 
constraints  imposed  by  the  known  projection  data.  The  reason  the  POCS  is  used  here  is  that,  even  in  the 
case  of  sparse  sampling,  the  size  of  the  projection  data  sets  can  be  large,  and  POCS  can  efficiently  handle 
large  data  sets. 

4  Results 

We  have  conducted  numerical  studies  to  investigate  image  reconstructions  by  using  the  TV-minimization 
algorithm  from  tomosynthesis  data  acquired  with  planar  and  non-planar  trajectories  shown  in  Fig.  1.  We 
use  ro(A)  to  denote  a  source  trajectory.  In  our  study,  we  first  consider  classical  tomosynthesis  in  which 
the  source  trajectory  is  a  segment  of  the  planar  circular  trajectory,  as  shown  in  Fig.  la.  In  the  study, 
we  have  chosen  R  =  7.0  cm,  and  A  G  [ — tt/6,  7t/6]  .  A  flat-panel  detector  is  placed  perpendicular  to  the 
line  connecting  the  source  and  the  rotation  center.  The  source-to-detector  distance  is  10.0  cm.  A  3D 
phantom  is  used  to  generate  cone-beam  projection  data.  It  should  be  pointed  out  that  we  have  generated 
analytically  the  data  so  that  the  data  contain  the  so-called  continuous-to-discrete  inconsistency  to  reflect 
realistic  imaging  conditions  in  practical  tomosynthesis.  Such  inconsistency  may  have  a  significant  impact 
on  reconstruction  accuracy.  In  top  row  of  Fig.  2,  we  display  the  phantom  images  within  planes  specified 
by  x  =  0  cm,  y  —  0  cm,  and  2  =  0.5  cm,  respectively.  We  first  generated  projection  data  from  the 
phantom  at  30  views  uniformly  distributed  on  the  planar  circular  trajectory  described  above.  From  the 
data,  we  subsequently  reconstructed  images  by  using  the  TV-minimization  algorithm,  and  we  displayed  in 
the  middle  row  of  Fig.  2  the  reconstructed  images  within  planes  specified  by  x  —  0  cm,  y  =  0  cm,  and 
2  =  0.5  cm. 

We  also  consider  a  non-planar  trajectory  that  consists  of  two  orthogonal,  circular  segments,  as  shown 
in  Fig.  lb.  In  this  study,  geometric  parameters  such  as  R  and  detector-to-source  distance  were  chosen  to 
be  identical  to  those  in  the  above  study  for  the  planar  circular  trajectory.  However,  the  angular  ranges 
of  the  two  orthogonal  circular  segments  are  specified  by  A  G  [— 7r/6,7r/6]  and  A  G  [7r/6,7r/2],  respectively. 
Using  this  non-planar  trajectory,  we  generate  analytically  projection  data  from  the  same  phantom  over 
total  30  views  of  which  15  views  are  uniformly  distributed  on  each  of  the  two  orthogonal  circular  segments. 
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The  number  of  views  over  one  of  the  two  trajectories  is  only  one  half  of  that  for  the  planar  trajectory. 
Therefore,  the  total  imaging  radiation  doses  in  the  studies  involving  the  planar  and  non-planar  trajectories 
are  the  same.  From  the  generated  data,  we  reconstructed  images  by  use  of  the  TV-minimization  algorithm. 
In  the  bottom  row  of  Fig.  2,  we  show  the  reconstructed  images  within  planes  specified  by  x  =  0  cm,  y  =  0 
cm,  and  2  =  0.5  cm,  respectively.  We  also  display  in  Fig.  3  the  profiles  in  the  reconstructed  images 
along  the  axes  specified  by  x  =  0  cm  and  y  —  0  cm,  by  y  =  0cm  and  z  =  0  cm,  and  by  x  =  0cm  and 
2  =  0cm,  respectively.  Comparison  of  the  results  in  Figs.  2  and  3  suggests  that  tomosynthesis  with  a 
non-planar  trajectory  may  yield  more  data  information,  and  thus  images  with  higher  quality,  than  classical 
tomosynthesis  with  a  planar  trajectory. 


5  Conclusion 

In  the  work,  we  have  investigated  non-planar  trajectories  in  tomosynthesis  for  yielding  more  data  infor¬ 
mation  than  classical  tomosynthesis  with  a  planar  trajectory.  Our  preliminary  study  indicates  that  with 
the  same  number  of  projection  views  (or,  equivalently,  the  same  dose  level,)  the  new  imaging  configuration 
may  lead  to  images  with  improved  quality  over  that  obtained  with  classical  tomosynthesis.  The  proposed 
imaging  strategy  and  TV-minimization  reconstruction  algorithm  may  find  applications  in  IGRT,  security 
scan,  industrial  imaging,  and  sample/specimen  evaluation. 
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Figure  2:  True  image  (top  row)  and  images  reconstructed  from  30- view  projection  data  acquired  with 
the  planar  trajectory  (middle  row)  and  the  non-planar  trajectory  (bottom  row).  The  first,  second,  and 
third  columns  represent  the  images  within  2D  slices  specified  by  x  =  0  cm,  y  —  0  cm,  and  £  =  0.5  cm, 
respectively.  The  display  gray  scale  is  [0,  2].  The  horizontal  and  vertical  axes  have  a  unit  of  cm. 
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Figure  3:  Profiles  of  the  images  reconstructed  from  30-view  projection  data  acquired  with  the  planar 
trajectory  (dashed  curve)  and  the  non-planar  trajectory  (solid  curve)  along  (a)  the  z-axis,  (b)  the  x-axis, 
and  (c)  y- axis,  respectively.  The  dotted  curve  indicates  the  true  profiles  along  these  lines. 
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ABSTRACT 

Current  dedicated,  cone-beam  breast  CT  scanners  generally  use  a  circular  scanning  configuration  largely  be¬ 
cause  it  is  relatively  easy  to  implement  mechanically.  It  is  also  well-known,  however,  that  a  circular  scanning 
configuration  produces  insufficient  cone-beam  data  for  reconstrucing  accurate  3D  breast  images.  Approximate 
algorithms,  such  as  FDK  has  been  widely  applied  to  reconstruct  images  from  circular  cone-beam  data.  In  the 
FDK  reconstruction,  it  is  possible  to  observe  artifacts  such  as  intensity  decay  for  locations  that  are  not  within  the 
plane  containing  the  circular  source  trajectory.  Such  artifacts  may  potentially  lead  to  false  positive  and/or  false 
negative  diagnosis  of  breast  cancer.  Non-circular  imaging  configurations  may  provide  data  sufficient  for  accurate 
image  reconstruction.  In  this  work,  we  implement,  investigate  innovative,  non-circular  scanning  configurations 
such  as  helical  and  saddle  configurations  for  data  acquisition  on  a  dedicated,  cone-beam  breast  CT  scanner,  and 
develop  novel  algorithms  to  reconstruct  accurate  3D  images  from  these  data.  A  dedicated,  cone-beam  breast  CT 
scanner  capable  of  performing  non-circular  scanning  configurations  was  used  in  this  research.  We  have  investi¬ 
gated  different  scanning  configurations,  including  helical  and  saddle  configurations.  A  Defrise  disk  phantom  and 
a  dead  mouse  were  scanned  by  use  of  these  configurations.  For  each  configuration,  cone-beam  data  were  acquired 
at  501  views  over  each  turn.  We  have  reconstructed  images  using  our  BPF  algorithm  from  data  acquired  with 
the  helical  scanning  configuration. 


1.  INTRODUCTION 

There  has  recently  been  a  renewed  interest  in  developing  dedicated  breast  CT  scanner  because  of  its  potential 
to  provide  3D  images  of  high  quality,  potentially  leading  to  improved  sensitivity  and  specificity  in  breast  cancer 
detection.  Current  dedicated,  cone-beam  breast  CT  scanners1,2  use  a  circular  scanning  configuration  for  data 
acquisition.  It  is  well-known,  however,  that  a  circular  scanning  configuration  cannot  yield  cone-beam  data 
sufficient  for  reconstructing  accurate  3D  images  of  the  breast.  We  hypothesis  that  advanced,  non-circular  scanning 
configurations  and  appropriate  reconstruction  algorithms  can  yield  images  with  quality  higher  than  that  obtained 
with  a  circular  scanning  configuration  and  the  approximate  FDK  algorithm.  In  this  work,  we  investigate  and 
develop  advanced,  non-circular  scanning  configurations  for  data  acquisition  in  dedicated  breast  CT,  and  we  also 
develop  algorithms  for  reconstructing  accurate  3D  images  from  the  acquired  data. 

2.  METHOD  AND  MATERIALS 

2.1  Imaging  configuration 

In  this  study,  we  have  used  a  dedicated,  cone-beam  breast  CT  scanner  capable  of  performing  non-circular 
scanning  configurations  at  UC-Davis.35  With  the  patient  lying  prone  on  a  table  of  the  scanner,  the  patient’s 
breast  hangs  in  a  pendulant  position  through  an  opening  in  the  table.  The  source  and  flat-panel  detector  can 
be  rotated  around  and/or  translated  with  respect  to  the  patient’s  breast  for  collecting  cone-beam  data.  Using 
a  Defrise  disk  phantom  and  a  dead  mouse  as  the  imaged  subjects,  we  have  investigated  several  advanced,  non¬ 
circular  scanning  configurations  such  as  the  helical  and  saddle  scanning  configurations.  The  helical  configuration 
has  most  widely  been  adopted  in  diagnostic  CT  for  achieving  rapid  data  acquisition.  Helical  scan  is  a  natural 
extension  of  circular  scan,  and  it  has  been  shown  to  yield  data  sufficient  for  accurate  3D  reconstruction.  The 
saddle  scanning  configuration  can  be  interpreted  as  a  generalization  of  the  circular  scanning  configuration,  in 
which  the  z-component  of  the  X-ray  source  will  now  vary  periodically.  It  can  yield  data  sufficient  for  accurate 
3D  reconstruction  of  breast  images.  The  implementation  of  a  saddle  configuration  requires  the  X-ray  source 
movement  along  the  z-axis  according  to  a  non-linear  function  of  the  rotation  angle. 
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We  use  the  Defrise  disk  phantom  in  our  study  because  it  provides  a  test  on  the  reconstruction  accuracy 
along  the  longitudinal  direction.  For  each  configuration,  we  acquired  cone-beam  data  at  501  views  uniformly 
distributed  over  a  turn.  In  certain  non-circular  scanning  configurations  under  study  such  as  the  saddle  configu¬ 
ration,  the  X-ray  source  needs  to  be  translated  as  a  non-linear  function  of  the  source  rotation  angle.  Because  of 
the  mechanical  limitation  of  motor  control,  such  a  non-linear  movement  was  realized  through  a  series  of  analyti¬ 
cally  pre-determined  source  positions  and  additional  interpolated  positions  between  two  adjacent  pre-determined 
source  positions.  In  the  implementation  of  the  saddle  scanning  configuration,  we  computed  250  source  positions 
distributed  equally  over  2i r  using  an  analytic  formula.  Based  upon  the  250  pre-calculated  source  positions,  we 
interpolated  additional  250  source  positions  for  additional  data  collection. 

2.2  Reconstruction  algorithm 

In  the  last  few  years,  significant  advances  have  been  made  in  the  development  of  theoretically  exact  algorithms 
for  image  reconstruction  from  cone-beam  projections.  These  algorithms  were  initially  developed  for  a  helical 
trajectory.  However,  they  have  subsequently  been  extended  to  reconstruct  images  from  cone-beam  projections 
acquired  with  other  trajectories,  for  example,  saddle  trajectory  ,  line+circle  trajectory  and  other  general  tra¬ 
jectory.  One  of  these  algorithms  is  the  back-projection  filtration  (BPF)  algorithm  for  reconstructing  images  on 
chords.6  A  unique  property  of  the  BPF  algorithm  is  that  it  can  reconstruct  an  image  on  a  chord  from  truncated 
data  as  long  as  the  data  on  the  intersection  segment  of  the  chord  with  the  object  is  non-truncated.  It  has  been 
shown  that  the  property  of  the  BPF  algorithm  can  be  exploited  for  obtaining  an  region  of  interest  (ROI)  image 
from  projection  data  containing  truncations.  ROI  imaging  can  reduce  the  dose  to  patient  and  it  has  the  potential 
benefits  of  scatter  reduction.  The  reconstruction  formula  for  BPF  can  be  written  as: 


j  1  P3'c2  _ 

fc(xc,Xa,Xb )  =  x~2  n -  u  =[  - c—\f(xc  2  -  X'c)(x'c  -  x^jgc(x'c,  Xa,  A;,)  +  C]  (1) 
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where  A  denotes  the  rotation  angles  of  the  x-ray  source,  chord  is  defined  as  the  line  segment  connecting  two 
points  ro(Aa)  and  ro(A^)  on  the  trajectory,  xc  denotes  a  point  on  the  chord,  which  relates  to  the  spatial  vector  as 
r(xc)  =  ^[rb(Aa)  +  ro(Ab)]  +xcec,  ec  indicates  the  direction  along  the  chord,  xc  G  [ — Z ,  Z] ,  Z  =  \  |rb(Aa)  —  ro  (A&)| 
denotes  one  half  of  the  chord  length,  /c(xc,  Aa,  Xb)  is  the  image  function  on  the  chord,  gc(xc ,  Aa,  A&)  is  the  back- 
projection  of  weighted  data  derivative  onto  the  chord,  (3  is  the  unit  vector  along  the  projection  direction,  and 
ew  points  along  the  rotation  axis.  The  detector  is  assumed  to  be  a  flat  panel  detector,  and  P(ud,Vd,  X)  is  the 
projection  data  at  view  angle  A. 


The  helical  trajectory  under  study  can  be  expressed  as 


f0(A) 


R  cos  A,  R  sin  A,  —  A 
27 r 


(3) 


where  R  is  the  distance  from  the  source  point  to  the  rotation  axis,  and  h  is  the  pitch  of  the  helical  trajectory, 
which  is  defined  as  the  translation  distance  of  the  imaged  object  during  one  turn  of  gantry  rotation. 

The  saddle  trajectory  under  study  can  be  expressed  as 

rb  (A)  =  ( R  cos  A,  R  sin  A,  h  cos  2A)T  (4) 

which  can  be  interpreted  as  a  generalization  of  circular  trajectory. 
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(a)  (b)  (c) 

Figure  1.  (a)  Image  reconstructed  by  use  of  the  FDK  algorithm  from  circular  cone-beam  data,  (b)  Image  reconstructed  by 
use  of  the  BPF  algorithm  from  helical  cone-beam  data,  (c)  Image  reconstructed  by  use  of  the  BPF  from  saddle  cone-beam 
data.) 


(a)  (b)  (c) 

Figure  2.  Images  within  (a)  transverse,  (b)  sagittal,  and  (c)  coronal  slices,  respectively,  reconstructed  by  use  of  the  BPF 
algorithm  from  real  data  acquired  with  a  helical  scanning  configuration. 


3.  RESULTS 

We  have  conducted  studies  using  both  computer-simulated  and  real  phantom  data.  In  the  simulation  studies, 
we  acquire  cone-beam  data  at  501  views  uniformly  distributed  over  2tt  angular  range.  A  flat-panel  detector  used 
has  a  size  of  512x512,  and  an  image  is  reconstructed  on  a  256x256x256  array.  Using  the  Defrise  phantom,  we 
first  generated  cone-beam  projection  data  at  501  views  with  the  circular  trajectory.  From  the  circular  cone-beam 
data,  we  reconstructed  images  by  using  the  FDK  algorithm.  In  Fig.  la,  we  display  one  of  such  reconstructed 
images  in  which  strong  artifacts  can  be  observed.  We  also  generated  cone-beam  projection  data  using  a  helical 
trajectory  and  a  saddle  trajectory.  From  these  data,  we  used  the  BPF  algorithm  to  reconstruct  images,  which 
are  shown  in  Figs,  lb  and  lc. 

We  have  also  implemented  helical,  saddle,  and  other  non-circular  scanning  configurations  on  the  dedicated 
breast  CT  scanner.  Using  these  configurations,  data  were  collected  from  the  Defrise  disk  phantom  from  which 
we  reconstructed  images  by  use  of  the  BPF  algorithm.  In  Fig.  2,  we  show  the  reconstructed  images  within  (a) 
the  transverse,  (b)  sagittal,  and  (c)  coronal  slices,  respectively,  from  data  acquired  with  the  helical  scanning 
configuration.  Clearly,  the  helical  scanning  configuration  can  yield  data  sufficient  for  accurate  reconstruction  of 
3D  images.  Moreover,  we  used  the  helical  configuration  to  acquire  a  data  set  from  the  mouse  from  which  we 
reconstructed  a  3D  mouse  image  using  the  BPF  algorithm.  In  Fig.  3,  we  display  mouse  images  within  (a)  a 
transverse,  (b)  a  sagittal,  and  (c)  a  coronal  slice,  respectively. 
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(a)  (b)  (c) 

Figure  3.  Mouse  images  within  (a)  transverse,  (b)  sagittal,  and  (c)  coronal  slices,  respectively,  reconstructed  by  use  of  the 
BPF  algorithm  from  data  acquired  with  a  helical  scanning  configuration. 


4.  CONCLUSION 

In  this  work,  we  have  conducted  a  preliminary  investigation  and  implementation  of  non-circular  scanning  con¬ 
figurations  on  a  dedicated  breast  CT  scanner  for  acquiring  cone-beam  data  sufficient  for  accurate  image  recon¬ 
struction.  We  also  studied  image  reconstruction  from  data  acquired  with  the  non-circular  trajectories  under 
studies  by  using  the  recently  developed  BPF  algorithm.  Results  of  the  studies  suggest  that  non-circular  scanning 
configurations  can  be  developed  on  dedicated,  cone-beam  breast  CT  for  acquiring  data  sufficient  for  accurate  3D 
reconstruction  of  breast  images. 
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